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FOREWORD 


This  formal  technical  report  entitled  "Ground  Motion 
Predictive  Techniques  for  Porous  Geologic  Media,"  is  sub- 
mitted by  Systems,  Science  and  Software  (S3)  to  the  Advanced 
Research  Projects  Agency  (ARPA)  and  to  the  Defense  Nuclear 
Agency  (DNA) . The  report  presents  the  results  of  the  third 
phase  of  a continuing  effort  to  develop  reliable  material 
models  and  computer  techniques  for  predicting  the  motion  of 
inhomogeneous  and  porous  geologic  media.  This  work,  in 
support  of  the  PRIME  ARGUS  and  MILITARY  GEOPHYSICS  programs, 
was  accomplished  under  Contract  No.  DASA  01-69-C-0l59(P00003) , 
which  was  funded  by  ARPA  and  monitored  by  DNA.  Dr.  Stanley 
Ruby  was  the  ARPA  Program  Manager  and  Mr.  Clifton  B.  McFarland 
was  the  DNA  Project  Scientist. 

Dr.  T.  David  Riney  was  the  S3  Project  Manager  for  the 
study.  The  technical  results  presented  in  this  report  repre- 
sent the  work  of  a number  of  S3  staff  members  in  addition  to 
the  authors.  It  is  appropriate  to  list  here  the  contributors 
to  technical  Sections  II  through  V. 


Section 

II : 

J. 

W. 

Kirsch,  A.  J.  Good,  G.  D.  Anderson 

Section 

III : 

J. 

K. 

Dienes,  G.  D.  Anderson, 

K. 

G. 

Hamilton 

Section 

IV: 

S. 

K. 

Garg,  D.  H.  Brownell, 

R. 

J. 

Archuleta 

Section 

V: 

G. 

A. 

Frazier,  S.  K.  Garg 

Dr.  M.  H.  Rice  participated  as  a consultant  on  this  project. 


1 


CONTENTS 


Foreword. 

Abstract. 


Page  No. 

. i 
• v 


I.  INTRODUCTION 


HOMOGENIZED  TREATMENT  OF  POROUS  WET  TUFF 


2.1  Introduction. 


2.2 

Pressure-Equilibrium  Mixtures  . 

/ 

1.1 

2.2.1 

Void-Free  Mixtures  . . . 

11 

2.2.2 

Mixtures  with  Voids.  . . 

14 

2.2.3 

S3  Mixture  Crush  Model 
(Disconnected  Pores)  . . 

19 

2.3 

Tabular  Arrays  of  Mixture  Equations 
of  State.  ..... 

25 

2.4 

Check 

Out  Calculations  with  TAMEOS 

29 

2.4.1 

Equations  of  State  of 
Constituents --Water  and  Tuff  . . . 

29 

2.4.2 

Check  Out  Calculati ons - - 
Planar  SKIPPER  .... 

31 

2.4.3 

Check  Out  Calculation-- 

Spherical  SKIPPER 

32 

2.5 

Parametric  Studies  of  Tuff/Water  Mixtures 

37 

2.5.1 

Source 

38 

2.5.2 

Peak  Radial  Stress  . . . 

38 

2.5.3 

Cavity  Growth 

41 

2.5.4 

Stress  vs_  Time  Profile  at  R = 40  m 

44 

2.5.5 

Ground  Motion  of  Media  at  R = 40  m 

46 

IMPROVED  PREDICTIVE  METHODS  FOR  GRANITE 

51 

3.1 

Introduction 

51 

3.2 

Finite 

Deformation  Theory  .... 

53 

L 


11 


Contents  (continued) 

Page  No. 

3.3  Equations  of  Motion 

3.4  Constitutive  Equations  for  Strain 

Hardening  Materials  58 

3.4.1  Equation  of  State 59 

3.4.2  Hardening  Models  for  Rock  Behavior  . 60 

3.4.3  Flow  Law  for  Materials  with 

Hardening 

3.4.4  Calculation  of  the  Multiplier.  ...  73 

3.4.5  Temperature  and  Strain  Rate  Effects.  77 

3.5  Set-Up  of  Spherical  Explosion  Calculation 

in  Granite 82 

3.6  Underground  Shot  Data 8 7 

3.7  Comparison  of  Calculations  and  Field  Data  . 90 

3.7.1  Cap  Model  Calculation 90 

3.7.2  Kinematic  Hardening 92 

3.7.3  Comparisons  of  Calculations 100 

3.8  Parametric  Studies  of  Granite  103 

3.9  Discussion  of  Results 122 

IV.  THEORY  OF  INTERACTING  CONTINUA  129 

4.1  Introduction 129 

4.2  General  Conservation  Laws 131 

4.3  Interaction  Terms 13A 

4.3.1  Interaction  Energy  Term:  p \p  ...  . 138 

4.3.2  Relation  to  Earlier  Work 141 

4.4  Constitutive  Relations 

4.4.1  Improved  Formulation  for 

Deviatoric  Stresses 142 

4.4.2  Thermodynamic  Constitutive  Laws.  . . 148 

4.4.3  Crushup  Model 149 


Contents  (continued) 


V. 


VI 


Page  No. 
. 152 


4.5  Thermodynamic  Porous  Code 

4.5.1  Conservation  Law  Equations 

m Moving  Coordinates 

4.5.2  Finite  Difference  Equations 

4.6  Stress  Pulse  Propagation.  . . 

* • • • 

FO^LCGGERrNGREARraS^fTI°N  AS  " MECHANISM 

* • • • 

5.1  The  Potential 


brAUerintrihef°?:rTriF-ering  Earth^akes 
7 Altering  the  Ground  Water  Conditions 

5.2  Alternate  Theoretical  Formulations. 


5.3 

Quasistatic  Linear  TINC 

Assumptions  . . 

5.4 

Conservation  Equations. 

5.5 

Constitutive  Equations. 

5.6 

Fluid  Seepage  

5.7 

Elastic  Deformations  in 

the  Rock  Matrix  . 

5.8 

FRI  Code.  . . 

5.9 

Injection  Well.  . 

5.10  Summary  and  Conclusions 
DISCUSSION  ... 


. 166 

. 173 

. 173 
. 176 
. 177 
. 178 
. 180 
. 188 
193 
195 
200 
209 

211 


REFERENCES. 

• • 215 

APPENDIX  A: 

SOUND  SPEED  IN  A POROUS  MATERIAL 

• • 223 

APPENDIX  B: 

H.E.  TEST  PARAMETER  CALCULATIONS 

• • 227 

APPENDIX  C: 

INTERACTION  TERMS. 

• • 241 

APPENDIX  D: 

EFFECTIVE  STRESS  LAWS 

• • 251 

IV 


ABSTRACT 


The  work  reported  consists  of  three  task  areas: 

(a)  development  of  constitutive  models  and  computer  methods 
for  calculating  stress  wave  effects  in  geologic  media  in  the 
vicinity  of  a buried  energy  source,  (b)  verification  of  the 
computer  models  and  their  application  to  examine  the  sensi- 
tivity of  ground  motion  predictions  to  the  material  parameters 
assumed  in  the  constitutive  models,  and  (c)  development  of 
methods  for  calculating  the  perturbation  of  residual  tectonic 
stress -strain  distributions  induced  by  changing  the  pore 
water  pressure.  A general  computer  subroutine  (TAMEOS)  is 
described  which  generates  thermodynamic  equations  of  state 
of  porous  wet  media  for  use  with  a table  look-up  procedure 
in  standard  ground  motion  codes.  TAMEOS  has  been  applied  to 
NTS  tuff  and  is  used  in  the  SKIPPER  code  for  a series  of 
spherical  calculations  in  which  the  crushup  strength  and 
volume  fractions  of  the  rock-water-void  mixture  were  varied. 
Four  constitutive  models  and  associated  subroutines  have 
been  incorporated  into  SKIPPER  for  improved  ground  motion  cal- 
culations for  rocks  with  high  shear  strength.  The  cap 
model  is  generalized  to  treat  the  full  range  of  pressure  and 
strain  encountered  in  underground  tests.  Generalized 
Mohr -Coulomb  models  include  one  without  work  hardening,  one 
with  isotropic  work  hardening,  and  one  with  kinematic  work 
hardening.  Model  sensitivity  calculations  for  granite  are 
presented  for  tne  four  SKIPPER  options.  Detailed  comparison 
of  calculations  with  field  measurements  are  presented  for 
the  kinematic  work  hardening  and  cap  models.  A thermodynamic 
formulation  of  the  Theory  of  Interacting  Continua  (TINC) 
is  presented  as  well  as  the  numerical  procedure  used  in  the 
new  POROUS  code  for  treating  spherical  ground  motion  problems 


v 
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in  the  TINC  framework.  Limited  calculations  using  the  new 
code  are  presented  for  partially  saturated  tuff.  Linearized 
TINC  equations  are  developed  for  describing  the  interaction 
of  a pore  fluid  with  a rock  matrix  as  the  fluid  is  driven 
through  the  geologic  mass  under  a hydraulic  gradient.  The 
2-D  finite  element  computer  code  fFRI)  for  treating  these 
geohydrological  processes  is  described.  Test  calculations 
for  the  rock-fluid  interactions  in  the  vicinity  of  a fluid 
injection  well  are  presented. 
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I.  INTRODUCTION 


Adequate  material  response  models  and  associated  com- 
putational techniques  are  required  if  ground  motion  predic- 
tions are  to  be  made  with  confidence.  One  is  concerned  with 
a signal  which  has  attenuated  from  stress  levels  in  the 
source  vicinity,  which  may  be  megabars,  to  stresses  at  large 
distances,  which  are  small  compared  to  strengths  of  earth 
media . 

In  this  transition  region  between  the  hydrodynamic 
source  and  the  distant  elastic  region,  the  material  response 
models  should  consider  such  complex  nonlinear  physical  pro- 
cesses as  dynamic  void  compaction,  heterogeneity,  pore  water 
pressure  and  diffusion,  yield  and  fracture  phenomena,  dilatancy, 
water  and  rock  interactions,  material  phase  changes,  and 
dependence  of  strength  parameters  on  the  thermodynamic  state. 
This  report  describes  improved  techniques  for  predicting  ground 
motion  that  have  been  developed  in  the  current  study.  The 
general  approach  followed  here,  as  in  earlier  work  described 
in  3SR-267  and  3SR-648,*  is  to  construct  material  models  of 
increasing  sophistication  from  available  material  properties 
data  and  to  develop  the  required  numerical  methods  to  evaluate 
stress  wave  phenomena  as  each  additional  effect  is  introduced 
into  the  model. 

In  addition  to  the  development  and  application  of  im- 
proved techniques  for  predicting  ground  motion,  the  report 
also  describes  work  directed  toward  understanding  earthquake 
triggering  mechanisms.  A computer  model  is  being  developed 
to  calculate  the  quasi -static  perturbation  in  a residual 
tectonic  stress-strain  distribution  as  a pressurized  fluid 
penetrates  through  a geologic  medium  under  a hydraulic 
gradient.  The  objective  is  to  analyze  the  effect  of  siting 

Here  and  throughout  this  document,  3SR-267  and  3SR-648  refer 
to  the  reports  [Refs.  1 and  2,  respectively]  describing  earlier 
phases  of  this  contract. 
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a reservoir  or  a deep  waste  disposal  well  in  the  vicinity  of 
a pre-existing  fault  zone. 

A typical  geologic  medium  consists  of  a rock  or  soil 
matrix  containing  cracks  or  pores  that  may  be  partially 
filled  with  water.  Even  if  the  matrix  material  is  unchanged, 
the  porosity  and  the  water  content  will  vary  with  depth  and 
with  surface  distance  and  the  stress  propagation  characteris- 
tics of  the  medium  will  vary  accordingly.  For  teleseismic 
calculations  it  is  impossible  to  know  the  porosities  and 
degrees  of  saturation  at  inaccessible  nuclear  test  sites. 

Even  when  local  geological  conditions  and  the  water  table 
location  have  been  established  by  field  logging  tests,  as 
would  be  possible  in  evaluating  the  vulnerability  of  under- 
ground structures,  it  is  economically  impractical  to  perforju 
laboratory  material  properties  tests  on  all  the  porosities 
and  degrees  of  saturation  that  occur.  Consequently,  it  is 
desirable  to  construct  the  material  models  in  such  a fashion 
that  the  response  of  the  medium  can  be  predicted  as  these 
quantities  are  varied. 

From  the  outset  of  the  present  study,  therefore,  the 
geologic  medium  has  been  considered  to  be  composite  and  a 
description  of  its  wave  propagation  characteristics  has  been 
sought  in  terms  of  the  behavior  of  the  isolated  rock  matrix 
and  water  components.  Reference  to  the  detailed  microstruc- 
ture of  the  composite  is  avoided,  however,  since  the  models 
are  to  be  used  in  continuum  type  computer  codes  in  which  the 
phenomena  of  interest  are  on  a much  larger  geometrical  scale. 
Consequently,  the  Theory  of  Interacting  Continua  (TINC)  was 
adopted  to  provide  a framework  general  enough  to  allow 
explicit  treatment  of  pore  pressure  effects  and  relative 
motion  between  the  rock  and  water  components.  Computer  codes 
usually  employed  for  ground  motion  calculations,  however, 
treat  a geologic  medium  as  a single  continuum  so  that  each 
incremental  volume  of  the  medium  has  associated  single 


values  of  pressure,  velocity,  etc.  Since  practical  calcula- 
tions are  currently  performed  using  such  codes,  material  response 
models  have  also  been  developed  to  be  compatible  with  applica- 
tion in  single  continuum  codes  by  adding  the  homogenizing 
assumption  of  no  relative  motion  between  the  rock  and  water. 

In  3SR-267  and  3SK-648  the  modeling  effort  centered  on 
a representative  tuff  at  the  Nevada  Test  Site  CATS  tuff)  to  be 
specific,  but  the  basic  methods  are  applicable  to  other  porous 
geologic  media  with  relatively  small  shear  strengths.  For 
example,  one  of  the  homogenized  composite  equation  of  state 
models  has  also  been  used  to  predict  the  behavior  of  clay 
shale  media  at;  the  MIDDLE  GUST  test  site.^  In  this  report 
the  model  development  has  also  specifically  treated  Cedar 
City  tonalite , a representative  geologic  media  with  a large 
shear  strength. 

Section  II  describes  an  homogenized  equation  of  state 
for  NTC  tuff  for  the  pressure  range  of  1 mbar  down  to  a few 
bars.  A computer  routine  has  been  developed  which  calculates 
the  isotropic  thermodynamic  states  of  rock-water- void  mixtu-es , 
including  a description  of  irreversible  collapse  of  the  air- 
filled  pores  (void  volume).  These  states  are  tabulated 
and  may  be  utilized  in  conjunction  with  a table  look-up  pro- 
cedure as  a subroutine  In  standard  ground  motion  codes. 

Primary  inputs  to  the  TAMEOS  subroutine  (for  Tabular  ^rrays 
Of  Mixtures  Equation  Of  State)  are  the  homogenized  model 
to  be  utilized  (e.g.,  one  of  the  PTEQ,  PEQ  or  P*EQ  models 
described  in  3SR-648),*  equations  of  state  of  the  isolated 
rock  and  water  components,  and  initial  volume  fractions  of 
rock,  water  and  air-filled  pores.  For  cases  in  which  experi- 
mental data  are  unavailable,  a simple  crushup  model  is 
employed  requiring  the  crushup  pressure,  sound  speed,  and 
elastic  crush  limit,  as  the  only  additional  inputs  to  be 


At  the  time  of  writing  TAMEOS  has  been  used  only  with  tables 
generated  with  the  PEQ  model,  but  may  be  readily  used  with 
other  pressure-equilibrium  models. 
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specified,  li  experimental  data  are  available,  the  crushup 
curve  can  be  directly  incorporated  into  the  TAMEOS  subroutine. 

The  TAMEOS  subroutine  has  be  )n  incorporated  into  the 
single  continuum  SKIPPER  code  and  a series  of  spherical  cal- 
culations made  using  the  PEQ  model  for  a representative  tuff 
with  varied  degrees  of  water  saturation  of  the  pore  space. 

An  8-kT  nuclear  source  is  simulated  by  a y = 1.4  gas  in  a 
cavity  of  initial  radius  of  3.72  m in  these  material  property 
sensitivity  calculation- . (In  Appendix  B the  results  of 
several  spherical  calculations  for  a 1000-lb  high  explosive 
energy  source  of  interest  to  the  NTS  underground  test  pro- 
gram are  also  presented.] 

Section  III  describes  the  basic  theory  and  the  asso- 
ciated series  of  subroutines  that  have  been  incorporated  into 
SKIPPER  for  Improved  ground  motion  calculations  for  high- 
shear-  strength  Igneous  rocks . ^Constitutive  laws  more 
general  than  the  Mohr-Coulomb  model  used  in  much  of  the  earlier 
work  are  introduced.  In  each  case,  an  associated  flow  rule 
is  assumed  in  developing  the  flow  law  and  the  influence  of 
finite  deformation  is  properly  treated  in  the  kinematic  rela- 
tions. Major  emphasis  was  on  the  work  to  generalize,  program 
and  test  the  Weidlinger  cap  model  (for  Cedar  City  tonalite ^ ) 
in  order  to  arrive  at  a constitutive  equation  that  covers 
the  full  range  of  pressure  and  strain  encountered  in  nuclear 
shots.  Generalized  Mohr-Coulomb  constitutive  equations  were 
also  incorporated  into  SKIPPER  as  options , including  one 
without  work  hardening,  one  with  isotropic  work  hardening, 
and  one  with  kinematic  work  hardening. 

These  four  modes  of  the  SKIPPER  code  were  exercised  in 
a series  of  material  model  sensitivity  calculations  for  a 
spherical  configuration  . These  parameter  studies  for  Cedar 
City  tonalite  are  all  for  the  same  simulated  8-kT  nuclear 
source  used  in  the  tuff  calculations.  Effects  of  the  strength 
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parameters  on  cavity  size,  stress  attenuation  and  displacement 
are  examined.  Most  of  the  work  described  in  Section  III, 
however,  concerns  comparison  of  the  Hardhat  and  Piledriver 
ground  motion  measurements  with  SKIPPER  calculations  using 
the  generalized  cap  model  and  the  kinematic  work  hardening 
model.  It  is  shown  that  the  latter  model  predicts  ground 
motion  in  much  better  agreement  with  field  measurements  than 
is  possible  when  the  cap  model  is  employed. 

A modification  of  the  TINC  formulation  to  more 
realistically  treat  the  deviatoric  stresses  in  the  rock 
matrix,  and  its  extension  to  account  for  thermodynamic 
effects,  is  presented  in  Section  IV.  The  1-D  computer  code 
POROUS  for  computing  stress  wave  effects  has  been  completely 
rewritten  to  account  for  the  more  comprehensive  TINC  model 
The  code  now  treats  spherical  as  well  as  planar  configurations 
w ereas  the  initial  version  for  solving  the  TINC  equations 
within  the  mechanical  formulation  treated  only  plane  waves 
(see  3SR-267  and  3SR-648] . Limited  calculations  using  this 
new  POROUS  code  for  partially  saturated  NT?  tuff  are 
presented. 


The  mechanical  Interaction  of  a pore  fluid  with  a 
saturated  rock  matrix  as  the  fluid  is  driven  through  the 
geologic  mass  surrounding  a fluid  injection  well  has  been 
modeled  within  the  linearized  TINC  equations.  The  formulation 
presented  m Section  V,  couples  the  deformation  and  diffusion 
fields  of  interest  to  the  problem.  In  order  to  solve  the 
associated  set  of  linearized  quasi-static  equations,  re- 
taining all  of  the  potentially  Important  interaction  terms 
a finite-element  method  of  solution  was  selected.  Available 
2-D  finite  element  computer  codes  for  treating  the  separate 
elastic  and  diffusive  processes  were  modified  and  combined 
into  a single  code  for  treating  the  coupled  processes.  This 
2-D  fluid-rock  interaction  code  (FRI)  is  described  in 
Section  V along  with  some  test  calculations. 
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In  Section  VI,  the  status  of  the  work  is  summarized 
and  suggestions  are  made  for  the  direction  of  the  effort 
during  the  next  contract  period. 

It  seems  appropriate  here  to  record  a number  of 
technical  publications  that  have  appeared  in  the  open  litera- 
ture describing  aspects  of  the  work  reported  in  3SR-267  and 
3SR-648.  Gurtman,  Kirsch  and  Hastings ^ presented  an  analy- 
tical equation  of  state  for  compressed  states  of  water. 

Morland  ^ described  the  initial  version  of  the  TINC  formula- 
tion  for  fluid  saturated  materials.  Garg L J presented  numeri- 
cal results  describing  wave  propagation  effects  using  the 
mechanical  TINC  formulation.  Garg  and  Kirsch^  showed  that 
the  TINC  framework  is  general  enough  to  include  various 
homogenized  composite  equations  of  state  (e.g.,  PTEQ,  PEO, 
and  P*EQ)  as  special  cases.  Norland^  presented  a finite 
deformation  plasticity'  theory'  with  isotropic  work  hardening. 
Additional  results  of  the  earlier  phases  of  the  work  have 
also  been  described  in  a number  of  oral  presentations  at 
technical  symposia. 
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II.  HOMOGENIZED  TREATMENT  OF  POROUS  WET  TUFF 


2.1  INTRODUCTION 


One  of  the  most  important  problem  areas  in  the  develop- 
ment of  ground  motion  codes  is  the  treatment  of  hydrodynamic 
rock/water  mixtures.  It  has  been  demonstrated  in  3SR-648 
that  one  may  derive  various  mixture  equations  of  state  on 
the  basis  of  a number  of  "equilibrium"  conditions  achieved 
behind  a shock  wave.  A unique  set  of  shock  states  is  speci- 
fied only  when  a constraint  is  prescribed  for  the  (shock) 
partitioning  of  internal  energy  between  the  rock  and  the 
wa t e r.  .Such  a set  of  states  can  be  obtained  if  the  pressure 
and  temperature  of  the  constituents  are  equal  (PTEQ) . If 
there  is  insufficient  time  for  thermal  equilibration,  but 
the  components  are  homogenized  to  the  point  that  they  are  in 
pressure  equilibrium,  other  shock  state  specifications  may  be 
made,  such  as  in  the  PEQ  and  P*EQ  models  discussed  in  3SR-648. 
The  latter  formulations  are  of  great  importance,  since  all 
laboratory  (and  mo.  t field)  experiments  designed  to  measure 
ground  motion  parameters  fall  within  the  regime  of  thermal 
nonequilibrium. 

Should  analytic  expressions  for  the  rock  and  water 
components'  equations  of  state  be  available,  it  is  clear 
from  the  preceding  study  (3SR-648)  that  an  analytic  for- 
mulation of  the  geologic  composite  does  not  result  for  the 
mixture  models  commonly 'utilized.  The  situation  is  further 
complicated  by  the  presence  of  voids  in  the  mixture.  Hence, 
it  should  be  recognized  that  the  simplification  associated 
with  modeling  the  effects  of  partial  water  saturation  by 
deriving  a single  hydrodynamic  equation  of  state  of  the  mix- 
ture, will  generally  require  tabular  arrays  of  mixture  states 
that  must  be  utilized  in  place  of  an  analytic  expression. 
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In  the  following  section  we  describe  a new  computa- 
tional tool,  called  Tabular  Arrays  of  Mixture  liquations  of 
otate,  IAMLOS.  Ibis  routine  calculates  homogenized  mixture 
states  and  stores  tnem  in  a tabular  array.  The  table  is 
stored  in  the  computer  and  individual  states  are  retrieved 
b>’  a rapid  table  look-up  routine.  In  its  simplest  form, 
t lie  table  consists  of  a rectangular  array  of  specific  volumes 
(V  J and  specific  internal  energies  (L)  and  tne  corresponding 
pressure  (p)  . To  treat  irreversible  pore  crushup  for  a 
partially  saturated  rock-water- void  mixture  a fourth  variable 
(a)  must  be  introduced  to  monitor  the  current  stage  of  the 
crushup  process.  In  the  w'ork  to  be  described  here,  the 
homogenized  treatment  of  the  porous  wet  mixtures  has  first 
been  specialized  to  the  PliQ  model  under  the  disconnected 
pores  postulate  fsee  3SR-648).  All  of  the  air-filled 
porosity  (void  space)  is  presumed  to  be  distributed  within 
the  rock  under  this  postulate.  The  porous  rock  and  water 
components  are  then  considered  to  be  in  pressure  equilibrium; 
the  two  components  ire  assumed  to  shock  to  the  same  states 
as  if  isolated  and  to  isentropically  release  witnout  any 
heat  transfer  between  them.  Ihis  formulation  was  selected 
ioi  first  treatment  because  of  its  relative  simplicity. 

(Other  pressure  equilibrium  formulations  could  be  tabulated 
in  a similar  manner.  It  is  planned  to  treat  the  P*LQ  and 
PTIiQ  models  described  in  5SR-048.) 

Tor  tnis  model  u is  defined  in  terms  of  tiie 
rock-water- void  volume  fractions  as  follows* 

(1)  (3) 

_ n + n 

a rrr~  (2>1) 

n 


Alternatively,  the  crushup  parameter  could  be  defined  in 
terms  of  the  distention  ratio  for  the  mixture 

(1)  (2)  (3) 

„ = n + n + n 

M (1)  (2) 

n + n 
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where 


U) 

n » rock  matrix  material  volume  fraction 
(3) 

n - void  volume  fraction  (air-filled  porosity)  . 

( (2)  (1)  (2)  (3)  \ 

\ n represents  the  water  volume  fraction,  n+n+n  =1/ 

Mixture  states  must  be  calculated  for  each  a > 1.  Hence, 
this  irreversible  crushup  regime  requires  a three-dimen- 
sional table  (V,  E,  and  a)  to  specify  p. 

As  in  the  case  of  a porous,  single  component  media, 
a suitable  expression  for  a(p,V)  is  required  to  utilize  the 
tabular  array  of  states  in  the  incompletely  crushed  regime. 

In  the  interest  of  generality,  the  crush  curve  may  be  speci- 
fically prescribed  if  enough  experimental  data  were  available. 
However,  for  calculational  purposes,  a specific  form  of 
crush  curve  for  the  rock  component  has  been  developed  for 
TAMEOS  which  incorporates  the  key  aspects  of  the  physics  of 
such  processes,  and  requires  a minimum  of  parameter  specifi- 
cations (crush  strength,  sound  speed  in  poreless  material, 
elastic  crush  regime  limit). 

The  mixture  states  have  been  computed  for  six  mix- 
tures of  tuff  and  water.  The  water  equation  of  state  in 
3SR-648  has  been  supplemented  to  include  the  high  pressure 
regime  (p  > 200  kbar  new  poreless  tuff  equation  of 

state  has  been  developed  especially  for  use  at  shock  pressures 
as  high  as  a megabar  SKIPPER  calculations  were  made  to  check 
out  the  TAMEOS  routines  and,  subsequently,  a parametric  study 
was  conducted  to  ascertain  the  effects  of  varying  water  and 
void  (air-filled)  volume  content  on  the  ground  motion  re- 
sulting from  8-kT  underground  nuclear  explosions.  The  volume 
fractions  of  tuff  matrix  material  considered  were  0.95  and 
0.8.  Three  degrees  of  water  saturation  were  considered  for 
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each;  all  pores  empty,  half  tne  pores  filled  with  water 
(.naif  saturated),  and  all  pores  filled  with  water  (fully 
saturated).  Quantitative  comparisons  of  tnese  computer 
calculations  are  presented  which  grapnically  illustrate 
the  effect  of  water  saturation  during  the  first  2U  msecs 
aitei  the  detonation.  Major  variations  were  observed  in 
the  stiess  time  iiistories  and  peak  stress  levels.  De- 
ereasing water  content  (and  therefore  higher  air-filled 
void  content)  leads  to  dramatically  lower  stress  levels, 
delayed  wave  arrival  times,  and  smaller  radial  displacements. 


Id 


L 


2.2  PRESSURE-EQUILIBRIUM  MIXTURES 


2.2.1  Void-Free  Mixtures 

Equations  of  state  utilized  in  ground  motion  codes 
usually  are  written  in  the  form 

p = p (V, E)  (2.2) 

Typically,  the  condition  of  pressure  equilibrium  between  the 
constituents  is  assumed,  i.e., 


where  the  i,j  subscripts  refer  to  the  ith  and  jth  component. 
Under  this  condition,  a hydrodynamic  pressure,  p » P.,  can  be 
assigned  to  the  mixture.  Hence,  overall  characterization  of 
the  mixture  is  obtained  by  taking 

p = P (V,E)  * Pi(v.,  Ej  = ...  (2.4) 

where 

k 

V * ^ M.Vi  (2.5) 

i = l 

k 

E = 2 MiEi  (2.6) 

i = l 

and  Mi  is  the  mass  fraction  of  the  ith  constituent. 

Obviously,  this  system  of  equations,  (2.4)  through  (2.6), 
is  not  determinate.  Given  a k component  mixture,  and  a set 
of  values  of  V and  E,  there  are  k+1  equations  for  2k 
unknowns  (V^  E.).  Thus,  k-1  additional  equati  ons  are  re- 
quired to  fully  specify  the  mixture  state. 


In  some  problems  of  interest,  there  may  be  enough 
time  for  the  constituents  to  thermally  equilibrate  (see  3SR-648 
for  a simplified  analysis  of  this  effect).  Under  this  condi- 
tion, one  can  close  the  algebraic  loop  by  requiring  the  tem- 
perature of  each  constituent  to  be  equal, 

T . = T . = . . . (2.7) 

1 j 1 J 

Of  course,  one  .nust  also  specify  the  caloric  equations  of 
state  of  the  constituents; 


The  pressure-thermal  equilibrium  blending  recipe  (PTEQ) 
given  by  Eqs.  (2.4)  through  (2.8)  can  be  solved  with  iterative 
computer  routines  for  the  simultaneous  solution  of  the  set 
of  algebraic  equations.  Generally,  these  take  considerable 
computer  time  and  in  hydrodynamic  calculations,  the  PTEQ 
formulation  can  most  efficiently  be  utilized  by  generating 
a tabular  array  of  V,  E states.  These  arc  used  in  place 
of  an  equation  of  state  with  a rapid  table  look-up  routine. 

Thermal  non- cqui 1 ibr ium  mixtures  arc,  conceptually,  more 
difficult  to  model.  The  PTEQ  formulation  can  be  treated  as  a 
homogeneous  material  whose  states  arc  uniquely  dclincd.  l\hen 
constituent  materials  arc  allowed  to  be  at  different  tempera- 
tures within  the  same  control  volume,  the  mixture  model  no 
longer  represents  a truly  homogeneous  material  (in  the  thermo- 
dynamic sense)  . 

The  PEQ  and  P*EQ  models,  introduced  in  3SR-648  and 
3SR-297,  arc  based  on  hypotheses  concerning  the  material 
interactions  under  shock  loading.  The  mixture  Ilugoniot  is 
derived  on  the  basis  of  these  assumptions.  Subsequent  re- 
lease states  arc  calculated  by  taking  pressure-equilibrium 
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mixtures  of  ti»e  release  adiabats  associated  with  the  shock 
states  of  each  constituent. 

In  the  PEQ  formulation,  each  constituent  is  assumed  to 
(shock)  compress  to  a state  on  the  Hugoniot  of  the  pure 
material.  Similarly,  upon  release  from  a shock  state,  each 
mixture  component  releases  along  the  isentropic  path  of  the 
pure  material.  Hence  the  shock  and  release  states  are  com- 
puted by  using  the  shock  pressure,  p^,  as  a parameter.  Each 
value  of  Pjj  implies  specific  values  of  , E^.  These  are 
then  utilized  in  Eq.  (2.5),  (2.6)  to  calculate  V,  E of  the 
PEQ  mix.  Since  these  values  satisfy  the  Hugoniot  relations, 
they  imply  that  certain  interactions  will  occur  which  result 
in  these  Hugoniot  states. 

Off-Hugoniot  states  are  obtained  by  calculating  the 
components'  release  isentropes  from  a given  shock  pressure 
and  imposing  pressure  equilibrium.  Thus,  the  interactions 
that  occurred  under  the  irreversible  shock-loading  process 
are  frozen  upon  release  and  each  constituent  expands  isen- 
tropically . 

One  may  take  a different  tack,  however,  and  specify 
that  interactions  occur  which  result  in  a mixture  Hugoniot 
that  differs  from  the  PEQ  version.  One  such  formulation  is 
the  P*EQ  model  introduced  in  3SR-648.  It  is  based  on  the 
assumption  that  the  entropy  of  the  water  component  is 
determined  by  a double-shock  process  and  remains  at  that 
level  behind  the  leading  portion  of  the  shock  wave.  The 
double-shock  results  from  the  impedance  mismatch  between 
the  rock  and  water. 

As  noted  in  3SR-648,  these  P*EQ  sets  of  water  and 
tuff  states  are  different  than  the  individual  Hugoniots  of 
each  constituent  that  are  used  in  the  PEQ  model.  This  is 


to  bo  expected  since  the  interaction  model  prescribes  lower 
values  of  water  entropy  at  a given  particle  velocity  than 
is  predicted  for  pure  water.  Once  these  sets  of  shock  states 
are  known,  isentropes  of  each  component  passing  through  a 
given  shock  pressure  may  be  put  into  pressure  equilibrium 
to  construct  the  mixture  release  adiubat  (as  in  the  Pi Q model). 
It  is  evident  that  the  PFQ  and  P*EQ  models  arc  very  similar 
in  concept.  They  differ  only  in  how  the  shock  interactions 
are  modeled.  Consequently,  the  mechanics  of  filling  a tabular 
array  of  states  is  exactly  the  same  for  either  mixture  model. 

2.2.2  Mixtures  with  Voids 

Accurate  predictions  of  the  response  of  naturally- 
occurring  geologic  materials  can  be  made  only  if  the  effects 
of  the  presence  of  air-filled  voids  (pores,  cracks,  etc.) 
are  taken  into  account.  In  the  case  of  shock  wave  propagation 
in  these  materials,  porosity  effects  are  manifest  in  two 
regimes:  pressure  levels  high  enough  so  that  material  is 
fully  crushed  and  the  lower  pressure  states  wherein  the 
voids  are  not  completely  removed.  In  both  regimes,  shock 
wave  propagation  is  retarded  due  to  the  extra  energy  required 
to  crush  the  material.  The  extra  shock  heating  due  to 
porosity  results  in  Ilugoniot  curves  which  arc  displaced  to 
the  right  of  the  poreless  material  Ilugoniot  in  the  p-V  plane 
(see  Fig.  2.1).  In  the  present  analysis,  the  mass  contri- 
bution due  to  the  compressed  air  has  been  ignored,  so  the 
excess  shock  heating  is  presumed  to  be  immediately  available 
to  the  rock/water  mixture. 

The  fully  crushed  regime  is  readily  adapted  to  the 
homogeneous  mixture  models  discussed  in  the  preceding  section. 
All  that  is  required  is  the  specification  of  a different  set 
of  constitutive  relations  to  partition  the  internal  energy 
behind  the  shock  wave.  In  PTEQ,  this  requirement  is  still 


Pressure  (kbar) 


0.28  0.32  0.36  0.40  0.44  0.48 

Volume  (cc/g) 


0.52 


0.56 


Fig.  2 . 1- -"Shifting"  of  the  Hugoniots  in  the  p-V 

plane  due  to  the  extra  shock  heating  associated 

with  initial  air-filled  porosity  for  two  PTEQ  NT'S 
tuff /water  mixtures.  (The  NT!>  tuff  equation  of 
state  is  given  in  3SR-648.) 


taken  care  of  by  the  assumption  of  thermal  equilibrium.  An 
alternative  model  is  the  initially  porous  version  of  the 
PEQ  model  (PEQP)  introduced  in  3SR-648.  It  maintains  that 
the  volume  fraction  of  voids  in  each  of  the  constituents 
is  equal  to  the  void  volume  fraction  of  the  mixture. 

Hugoniots  for  the  completely  crushed  constituents  can  be 
calculated  from  their  equations  of  state.  From  there  on, 
the  mixture  equation  of  state  calculational  procedure  is 
identical  to  that  outlined  for  the  PEQ  model.* 

The  energy  partition  is  directly  determined  by  the 
percentage  of  original  void  volume  that  is  assigned  to  the 
water  and  tuff  matrix  components.  Consequently , a multitude 
of  PEQ! * type  mixtures  exist,  each  based  on  different  sharing 
ol  the  void  volume  fraction.  The  disconnected  pores  postulate 
is  utilized  as  the  basic  mixture  model  in  this  report.  It 
assumes  that  all  voids  (air-filled  pores)  are  contained  in 
the  rock  matrix  and  that  the  bulk  pressure  in  the  distended 
tuff  is  in  equilibrium  with  the  effective  pressure  of  the 
water.  The  special  advantage  of  this  formulation  is  that  a 
minimum  of  experimental  information  is  required  to  charac- 
terize the  incompletely  crushed  regime. 

There  is  ample  evidence  that  porous  geologic  materials 
do  not  actually  lose  all  of  their  porosity  until  sufficiently 
high  compressions  are  achieved  [e.g.,  Refs.  IQ  through  12 ] • 
the  partially  saturated  rock  matrix  also  may  exhibit  precursor 
effects  which  do  not  appear  at  higher  degrees  of  water 
* ; 

The  derivation  of  mixture  equations  of  state  for  initially 
porous  materials  can  also  be  approached  from  different  van  - 
tage  points  such  as  that  used  in  developing  P*EQ  for  tuff/ 
water  mixtures.  In  this  instance,  the  water  is  presumed  to 
undergo  two  major  shock  waves  as  the  crush  wave  propagates 
through  the  material.  For  example,  the  leading  shock 
develops  from  the  initial  wave  and  the  secondary  shock  re- 
sults when  the  voids  are  closed.  In  3SR-648,  such  a model 
(P*EQP ) was  introduced  to  account  for  these  multiple  shock 
interactions. 
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saturation.  110  Hence,  to  account  for  this  behavior,  the 
crushing  of  the  partially  saturated  rock/water  composite 
should  include  an  elastic,  reversible  process  at  low  stress 
levels,  followed  by  a regime  wherein  plastic,  irreversible 
deformation  occurs  as  the  matrix  collapses.  These  are 
complex  processes  which  are  not  well  understood. 

A useful  model  has  been  suggested  for  metallic  foams 
by  Herrmann  that  simplifi  s discussion  of  the  problem. 

The  distension  parameter,  a,  is  assumed  to  be  a function  of 
stress  level, 

a = a(p)  (2.9) 

and  the  hydrodynamic  stress  is  computed  from  the  equation  of 
state  of  the  matrix  material,  with  the  modification  that  the 
effective  specific  volume,  V/ot,  be  utilized  in  place  of  V, 

i . e . , 

P = Pil  > E ) (2.10) 

Of  course, a(p)  is  not  generally  known.  It  is  only  a convenient 
functional  representation  that  Ignores  any  internal  energy 
effects.  Herrmann  discusses  some  general  forms  of  ot(p)  which 
have  been  employed  to  analyze  the  incomplete  crush  regime  of 
some  porous  metals.  Qualitatively,  these  are  sketched  in 
Fig.  2.2.  The  critical  parameters  in  such  a model  are  the 
pressure  at  which  all  pores  are  (irreversibly)  removed,  p , 
the  pressure  limit  to  the  elastic,  reversible  regime,  p , 
and  the  slope  of  a(p)  at  p = 0.* 

^ - 

a"(p)  at  p = 0 is  related  to  the  ratio  of  the  sound 
speed  in  the  porous  material  to  that  in  the  condensed  rock 
matrix  (see  Appendix  A) . 
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Fig.  2. 2- -Schematic  of  Herrmannls  model  for  the 
distension  ratio,  a(p) . Note  that  above  the 
elastic  limit,  release  from  a plastic  crush 
state  occurs  along  a prescribed  (reversible) 
release  path  which  could  result  in  some  void 
recovery  (a2  > a2).  Release  in  the  elastic 
regime  is  confined  to  the  loading  curve. 
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A review  of  the  tuff  crush-up  data  then  available 
was  reported  in  3SR-648.  A "universal"  crush  curve  for 
tuffs  was  derived  that  could  "fit"  the  available  crush  data 
when  adjusted  for  grain  density  variations  from  specimen  to 
specimen.  However,  the  error  band  in  this  universal  crush  law 
law  was  quite  large.  Hence,  if  experimental  crush  curves 
are  available  for  a given  tuff  (or  other  rock  matrix)  they 
should  certainly  be  utilized. 

The  S3  universal  fit  to  crush  data  was  functionally 
represented  by  a(p)  instead  of  Herrmann's  a(p)  formu- 
lation (where  p is  the  density).  The  concept  of  a crush 
density  was  easier  to  apply  to  the  data,  since  the  crush 
pressure  varied  considerably  for  the  different  tuffs.  However 
the  framework  created  by  invoking  the  disconnected  pores 
postulate  makes  it  especially  convenient  to  utilize  an  a(p) 
formulation  so  as  to  facilitate  the  derivation  of  pressure 
equilibrium  rock/water/void  mixtures.  Hence,  during  the  past 
year  a generalized  a(p)  crush  curve  formulation  has  been 
developed  for  the  rock  matrix  which  yields  similar  results  to 
the  universal  description  and  can  be  "tuned'^  to  fit  available 
experimental  information. 

2.2.3  S3  Mixture  Crush  Model  (Disconnected  Pores) 

The  basic  formulation  of  the  pressure  equilibrium  mix- 
ture with  porosity  in  one  component  (rock)  is  redefined  in 
terms  of  equating  the  bulk  pressure  of  the  porous  rock  matrix 
to  the  pressure  in  the  water.  As  pointed  out  in  3SR-648,  the 
bulk  pressure  in  an  isotropic  porous  rock  component  can  be 
defined  as 


where  the  subscript  1 denotes  the  poreless  rock  component. 
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That  is,  the  bulk  stress  exerted  on  a gage  is  predicted  to 
be  less  than  the  effective  pressure  assumed  in  Herrmann's 
formulation,  Eq.  (2.10).  Carroll  and  Holt[lsJ  recommend 
Eq,  (2.11)  in  preference  to  (2.10)  for  the  case  that  the 
pores  are  so  small  that  one  may  consider  the  rick  to  be 
homogeneous.  Another  argument  for  this  modification  can  be 
made  on  the  grounds  that  the  sound  speeds  predicted  by 

Eq.  (2.11)  are  physically  more  meaningful  (see  discussion 
in  Appendix  A) . 

The  algebraic  formulation  of  the  disconnected  pores 
mixture  model  is  given  by 

p * I T{sr  < 0 ' p2(v2.  e2)  • p(v,E) 

where 

V = M V + M V 
11  2 2 

E = M E +MV 

a - distention  ratio  of  rock,  some  function 
of  the  hydrostatic  pressure 

and  the  subscripts  1 and  2 denote  the  poreless  rock  and 
water  constituents  respectively. 

It  is  evident  from  Eq.  (2.11)  that  a prescription  of 
a is  all  that  is  required  to  completley  formulate  the 
porous  mixture  equation  of  state.  This  simplification  ig- 
nores any  modifications  to  a due  to  variations  in  moisture 
content  or  porosity.  (However,  a may  be  considered  as  an 
input  function  and  any  prescription  can  be  used  to  calculate 
the  states  of  a mixture.) 


(2.12) 

(2.13) 

(2.14) 

(2.15) 
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In  lieu  of  experimental  data  for  a particular  site,  a 
simple  set  of  expressions  have  been  developed  for  the  crushing 
of  the  porous  rock  constituent.  The  plastic  crushvregime 
has  been  modeled  by  a quadratic  expression: 

a = 1 * K-1))1  - '2-16> 

where 

ae  = distension  ratio  at  limit  of  elastic  region 

ie  = pressure  at  upper  limit  of  elastic  region 

Pc  = pressure  at  which  voids  are  completely  removed. 


The  quadratic  i 
and  da/dp  = 0 
a specification 
for  the  elastic 


s the  simplest  formulation  which  has  a = 1 
at  P = Pc  * Thus,  all  that  is  required  is 
of  the  crush  pressure  and  the  match  points 
region  (see  Fig.  2.3) . 


Release  from  any  point  in  this  region  occurs  with  cx 
fixed  to  the  minimum  value  achieved  during  loading.  Hence, 
to  derive  the  mixture  equation  of  state,  the  release 
states  are  computed  for  values  of  a in  the  vicinity 
of  the  Hugoniot  points  in  the  ct  > 1 regime  (see  discussion 
in  Section  2.3),  It  should  be  noted  that  a(p)  is  dependent 
on  the  direction  of  the  process,  and  is  not  a simple,  analytic 
function  of  pressure.  This  creates  the  need  for  a three- 
dimensional  table  in  this  regime,  so  that  with  a specified, 
valid  mixture  states  can  be  located  without  confusion.  The 
possible  overlap  of  release  states  in  the  p-V  plane  is 
readily  taken  into  account  by  this  method. 
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(kbar) 


Fig.  2 . 3- -Example  of  p-a  crush  curve  for  a sample  of 
51  porosify.  All  pores  are  completely  closed  at  1 kbar. 
Note  the  smooth  transition  across  the  elastic  limit 
(pe  - 0.2  kbar).  The  dashed  curve  represents  the  crush 
curve  without  an  elastic _ regime  (Pe  = 0).  The  initial 
siope  in  the  elastic  regime  was  chosen  to  produce  an 
initial  bulk  modulus  of  20  kbar  for  the  porous  mixture 
(compared  to  83.3  kbar  for  the  rock  material  and  8.9 
kbar  for  the  completely  plastic  curve. 
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Elastic  regimes  may  not  exist  for  all  mixtures  and 

plastic  crushing  may  begin  at  p » 0.  In  this  case,  p = 0, 

6 


(the  initial  porosity),  and  Eq.  (2.16)  reduces  to 


K -i)(J 


- E_ 

P. 


(2.17) 


However,  stress  wave  measurements  in  porous  samples 
have  indicated  the  presence  of  some  sort  of  precursor  phenomena. 
It  is  not  clear  that  one  can  characterize  this  as  an  elastic 
response  or  a str-.in  rate  effect.  Recognizing  the  potential 
need  for  a model  of  the  elastic  regime,  to  be  evaluated  in 
other  parametric  studies,  a formulation  has  been  included  in 
the  TAMEOS  routine. 

The  elastic  regime  in  the  S3  model  is  limited  by  a peak 
elastic  stress  level  (as  shown  in  Fig.  2.3).  Below  this  value, 
loading  and  unloading  takes  place  along  the  same  a(p)  curve. 

An  exponential  expression  for  ct(p)  in  the  elastic 

regime  was  selected  because  it  afforded  the  widest  choice 

of  parameters  to  match  to  sound  speed  and  bulk  modulus  data 

while  retaining  the  feature  that  \ . Tne 

V^P/P-O  \3p/P“Pe 

distention  ratio  is  given  by 


a * a ♦ 
o 


(a  ' ae)  I 

n-l  ) ' 


1 - e 


n P/P, 


(2.18) 


- n< 


(2.19) 


(°  ~ ae)  (P 

'P 


(k->) 


(2.20) 
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This  expression  is  smoothly  patched  (continuous  first  deri- 
vative) to  the  plastic  form  for  u(p)  at  p = p^. 

The  transcendental  relation,  Eq.  (2.19),  determines 
the  acceptable  sets  of  values  for  ag  and  pg.  Should 
sound  speed  and  bulk  modulus  data  be  available  for  the  porous 
rock  matrix,  an  additional  constraint  is  put  on  the  values 
of  a and  p . As  shown  in  Appendix  A,  the  two  elastic 

C w 

crush  parameters,  a and  p , are  related  to  the  sound  speed 

G G 

in  the  matrix  at  p = 0, 


c 


2 


Q 


f (da/ dp)  q ' 

1 - - K 

0 


(2.21) 


where  K , is  the  bulk  modulus,  c is  the  sound  speed  of 
0 0 

the  poreless  material,  and  (da/dp)Q  is  the  slope  of  the 

crush  curve  at  zero  pressure.  This  slope  determines  the 

degree  to  which  sound  waves  are  slowed  in  porous  materials. 

In  a stiff  matrix,  da/dp  * 0,  c = c . The  relationship  be- 

o 

tween  the  sound  speed  in  the  porous  matrix,  and  the  elastic 
parameters  is  obtained  by  differentiating  Eq.  (2.18)  and  sub- 
stituting into  (2.21)  to  arrive  at 


c = 


1 + 


c 4 
0 

HZ' 

pc-pe 


TT 


n 


(2.22) 
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2.3  TABULAR  ARRAYS  OF  MIXTURE  EQUATIONS  OF  STATE 


A computer  routine  has  been  written  to  generate 
tabular  arrays  of  thermodynamic  states  of  rock-water  mixtures. 
The  flow  chart  for  this  numerical  process  is  given  in  Fig.  2.4. 
Briefly,  let  us  consider  first  how  a table  is  constructed  for 
a saturated  rock-water  mixture.  For  given  equations  of  state 
of  the  two  materials,  the  PEQ  mixture  states  are  determined 
by  the  volume  fractions  of  each  constituent.  Additionally, 
a set  of  pressure  points  in  the  range  of  interest  and  the 
desired  mesh  fineness  for  V and  E must  be  specified. 

The  table  generating  program  then  computes  the  mixture  re- 
lease isentropes  from  states  on  the  Hugoniot  curve.  Based 
on  these  isentropes,  the  program  calculates  p-V-E  states  for 
values  of  V and  E which  are  suitable  for  the  table  look- 
up routine  incorporated  in  TAMEOS.  The  resulting  set  of 
pressure  points,  along  with  the  information  on  the  V and 
E meshes,  comprises  the  table  for  the  specified  (non-porous) 
mixture  in  the  given  range  of  pressures.  Other  states  in 
this  region  are  then  computed  in  TAMEOS  by  interpolation. 

If  ga_ -filled  pores  are  present  in  a given  mixture 
(a  = >1),  then  a separate  p-V-E  table  is  constructed  for 

each  a at  specified  mesh  intervals  in  the  range  [1,  a ] 

For  each  a there  is  a unique  pressure  point  on  the  porous 
rock  Hugoniot  given  by  the  solution  to 


where 

p = P (v  , E ) , equation  of  (2.24) 

1 1 state  for 

poreless  rock 
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Fig.  2. 4- -Flow  diagram  of  table  generating  portion  of  TAMEOS. 


and 


a = o(p)  , p < 


« = xt  P 1 Pc 


(2.25) 


There  is  only  one  isentrope  through  the  Hugoniot  at  this 
point,  and  to  allow  room  for  interpolation  in  the  table,  other 
isentropes  are  chosen  by  slightly  shifting  the  specific  volume 
to  either  side  of  the  Hugoniot  point  and  using  the  same  pressure. 
A table  of  p-V-E  states  is  then  obtained  for  each  a.  To  ob- 
tain the  pressure  at  a given  state  a-V-E,  the  two  tables 
corresponding  to  the  bracketing  values  of  a are  used  to  com- 
pute a pair  of  pressures  which  are  then  used  to  interpolate 
for  p at  the  given  a. 


In  most  situations,  the  desired  pressure  range  encom- 
passes both  porous  and  compacted  material  states.  Separate 
tables  must  then  be  constructed  for  the  two  regions,  the 
value  of  a determines  which  of  the  tables  to  use. 


Table  Look-Up 

The  pressures  which  comprise  a table  correspond  to 

values  of  V,  E,  and  a chosen  so  that  no  time  consuming 

search  is  needed  when  performing  an  interpolation  at  a given 

state.  Briefly,  the  scheme  consists  of  picking  values  for 

the  independent  variables  (V,  E,  and  a)  of  the  form 
k . k n 

2 + j*2  /2  , where  k may  be  any  integer,  n any  non-negative 

integer,  and  0 £ j < 2n.  By  varying  n one  can  control  the 
spacing  of  the  grid;  k and  j determine  not  only  the  values 
but  also  the  indexing  of  the  independent  variables.  Full 
advantage  is  then  taken  of  the  binary  representation  of  num- 
bers in  the  computer  to  expedite  the  index  calculation  as  well 
as  the  actual  interpolation  for  p . The  result  is  an  extremely 
fast  interpolation  with  a minimal  amount  of  arithmetic. 
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Application  of  TAMEOS  for  Porous  Mixtures 


In  calculations  employing  TAMEOS  a simplification  has 
been  introduced  for  the  crush  curve,  a(p) . Experience  has 
demonstrated  that  to  avoid  lengthy,  time  consuming  iterations 
for  a in  the  implicit  equation,  a = P | , E^jj, 

one  mav  reformulate  a(p)  as  a (fitted)  function  of  V . , 
along  the  Hugoniot,  i.e., 


[a(p)  ]„  = a (pH(vH))  = aR(Vmix) 


(2.26) 


This  reformulation  actually  simplifies  the  calculation  a 
great  deal;  given  a state  (Vm^x , Em^x)  , the  value  of  a 
can  be  quickly  determined  and  TAMEOS  furnishes  the  (inter- 
polated) pressure  for  pmix  at  (a,  Vmix , Enix) . Of  course, 
checks  are  required  to  deteuuine  whether  a cell  is  being 
loaded  or  undergoing  a release  process  so  that  a may  be 
computed  correctly. 
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2.4  CHECK  OUT  CALCULATIONS  WITH  TAMEOS 


A variety  of  calculations  were  conducted  to  verify  the 
projected  three-place  accuracy  of  the  TAMEOS  routine  and  to 
gain  experience  in  its  application  to  typical  ground  motion 
calculations.  In  this  regard,  the  retrieval  time  for  a value 
f°r  P*  given  E and  V,  was  determined  to  be  comparable  to 
that  of  typical  analytic  expressions  for  equations  of  state. 

It  is  difficult  to  quantify  this  difference  but  it  appears 
that  TAMEOS  can  be  utilized  as  quickly  as  any  single  component 
media  equation  of  state.  Of  course,  storage  of  the  tabular 
array  is  required. 

2.4.1  Equations  of  State  of  Constituents- -Water  and  Tuff 

The  water  equation  of  state  utilized  to  derive  rock/ 
water  mixtures  equations  of  state  is  composed  of  the  original, 
analytic  expressions  reported  in  3SR-648  and  a tabular  represen- 
tation for  states  with  entropies  greater  than  3.3*107  ergs/g  °K 
(corresponding  to  shock  pressures  higher  than  200  kbar).  The 
tabular  portion,  compiled  by  Bjork,^16^  is  valid  up  t'  10  Mbar 
in  the  vicinity  of  the  Hugoniot  and  treats  expansions  into 
the  gas  phase.  In  the  present  calculations,  the  bulk  of 
the  mixture  states  were  of  peak  pressures  below  200  kbar,  so 
that  the  water  equation  of  state  in  3SR-648  describes  most  of 
the  water  states  involved  in  this  series. 

A new  equation  of  state  for  poreless  tuff  has  been 
derived  which  is  based  on  the  Hugoniot  curves  for  saturated 
Rainier  Mesa  tuff  reported  by  Shipman,  et  al. ^ Points 
on  the  dry,  poreless  tuff  Hugoniot  were  obtained  from  this 
data  hy  assuming  that  the  tuff  is  in  pressure  equilibrium 
with  the  water  and  that  the  water  states  are  the  same  as 
those  postulated  for  the  pure  water  constituent.  (Hence, 
these  tuff  states  are  those  which  satisfy  the  PEQ  mixture 
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model,)  \t  each  pressure, 


VHM(p)  = MT  VHT(p)  + MW  VHW(p)  (2.27) 

where  VjM(p)  is  the  specific  volume  of  the  saturated  mixture 
at  pressure  p on  the  Hugoniot.  VHT(p)  and  V,IW(p)  are  the 

specific  volumes  of  the  tuff  and  water  constituents  at  the  same 
shock  pressure. 

It  was  determined  that  the  RTS  tuff  representation 
used  in  3SR-648  was  not  adequate  ft  these  higher  shock 
pressures.  A series  of  shock  velocity-particle  velocity  points 
were  computed  on  the  basis  of  new  tuff  Hugoniot  points  ex- 
tracted from  high  pressure  data.  These  values  did  not  lie  on 

a straight  line  but  could  be  fit  within  the  experimental 
accuracy  to  the  form 


U = 

where  ll  and  u are  the 

This  expression  implies  a 
by  [ 18,  19] 

PH(V)  = (pfln  r/2u  )£l  - (l 

P„(V)  = ( pfln  T/2u  )|l  + ( l 

where  n = (l  - V/V  \ and 

o ' 

u = d 2 7f4 


(2.28) 

velocities . 
Hugoniot  given 


- - \1/2  1 

4ya2 / A 2 | 

nb  < 1 

(2.29) 

U/21 

4ua2/p) 

nb  > l 

(2.30) 

(2.31) 

a + bu  + du2 

shock  and  particle 
pressure-densi ty 


A = (1  -bn)  2 - 2adrf2  (2.32) 
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The  two  branches  of  the  Hugoniot  curve  are  discussed  fully 
in  Ref.  18  and  19.  The  coefficients  in  Eq.  (7.28)  for  dr> 
compacted  Rainier  Mesa  tuff  are  a = 3.50  mm/psec,  b = 0.7047, 
and  c = 0.1005  (mm/psec)"1.  Since  in  compression  btT  < 1 
the  Hugoniot  for  this  tuff  will  always  be  represented  by 
Eq.  (2.27) 

The  Hugoniot  curves  are  used  as  the  basis  of  a p-V-E  equa- 
tion of  state  for  poreless  tuff  similar  to  that  adapted  by 
Butkovich, ^ ^ i.e.,  a Mie -Gruneisen  equation  of  state  with  the 

Gruneisen  ratio  proportional  to  specific  volume.  This  is  written 

P ' G„P,E  * PH-V’  l1  - (v,  - v)J  (2.33) 

where  G p is  the  product  of  the  Gruneisen  ratio  and 

Q 0 

density  at  normal  conditions.  It  was  determined  that 

Gopo  = 0,732  g/cm3  (3.34) 

provided  the  best  fit  to  release  data.  The  solid,  grain 
density,  , was  assumed  to  be  2.22  g/cm3. 

2.4.2  Check-out  Calculations - -Planar  SKIPPER 

As  a test  case,  the  initial  volume  fractions  were  se- 
lected to  simulate  a partially  saturated  tuff/water  mixture; 


(1) 
n = 

0 

O 

o 

tuff 

(2) 
n = 

0 

0.23, 

water 

(3) 

n = 
0 

O 

• 

o 

void . 

t 
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The  initial  distension  ratio  in  the  tuff,  a,  is  1.1  and  crush 
pressure  was  set  at  pc  = 5 kbar.  It  was  assumed  that  the 
mixture  had  no  strength  and  an  elastic  crush  regime  was  not 
included  ae  = a0 * Pe  " 0,  so  that  a is  given  by  Eq.  (2.17) 

Figure  2.5  is  a plot  of  the  pressure-volume  states 
calculated  with  the  planar  SKIPPER  code  for  a step-pulse  of 
Pr  = 4 kbar  followed  by  a release  wave.  The  curves  in  Fig.  2.5 
represent  theoretical  PEQ  Hugoniot  and  release  curves.  For 
the  crushup  wave,  the  computed  p-V-states  are  indicated  by 
the  squares,  and  the  release  states  are  represented  by 
circles.  The  loading  p-V  states  calculated  in  SKIPPER  are 
seen  to  lie  close  to  the  PEQ  Hugoniot.  Release  states 
follow  an  a = constant  locus  of  states  as  evidenced  by 
the  calculated  points  for  the  4-kbar  release  wave. 

Figure  2.6  shows  a similar  calculation  for  a 145-kbar 
step  pulse  (complete  crushing).  It  is  evident  that  the  low 
pressure  (p  < pc)  crushing  occurs  very  close  to  the  PEQ 
Hugoniot  p-V  trace,  while  the  higher  pressure  loading  p-V 
states  lie  to  the  right  of  the  Hugoniot  as  a consequence 
of  the  extra  internal  energy  contributed  by  the  q-term.  The 
q- term  insures  that  loading  is  along  the  Rayleigh  line  in 
order  that  the  shock  heat  at  p^  is  correct.  The  computed 
release  states  are  indeed  seen  to  be  in  excellent  agreement 
with  the  theoretical  PEQ  release  adiabat. 


2*4.3  Check-out  Calculations - -Spherical  SKIPPER 


The  parametric  study  to  investigate  porosity  effects 
on  wave  propagation  in  tuff  was  conducted  with  the  spherical 
SKIPPER  code.  These  calculations  are  discussed  in  detail 
in  the  next  section.  As  a check  on  the  accuracy  of  TAMEOS, 
(3) able  for  20*°  P°rous  dry  tuff  ( n = 0.80,  n = 0.00, 
n = 0 . 20 ) was  generated  and  utilized  in  an  identical  cal- 
culation to  that  made  with  the  analytical  equation  of  state. 
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Fig.  2 . 5- -Pressure-volume  states  taken  from  a (planar) 
SKIPPER  code  check-out  calculation  using  the  TAMEOS 
routine  in  the  plastic  crush  regime.  In  this  example 
a tuff/water  mixture  table  was  generated  and  used  in 
place  of  an  analytic  equation  of  state. 
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Fig.  2 . 6- -Pressure-volume  states  taken  from  a 
(planar)  SKIPPER  code  check-out  calculation  using 
TAMEOS  in  the  completely  crushed  region  (p  > p ) . 
The  shock  state  computed  with  TAMEOS  was  accurate 
to  three  places  in  this  example. 


(The  crushup  parameters  were  also  the  same.)  Figure  2.7  is 
a comparison  of  the  stress  wave  profile  computed  with  the 
analytic  and  tabular  equations  of  state  after  500  cycles. 
The  two  results  are  almost  identical  and  well  within  the 
projected  three-place  accuracy  of  the  table. 
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Fig.  2. 7- -Comparison  of  the  radial  stress  profiles  calculated 
with  the  analytic  tuff  equation  of  state  (Eq.,  2.33)  and  the 
tabular  array  of  states  for  a 20%  porous  tuff  media  with  no 
water  present  in  the  pores.  The  time  is  3.45  msec  after 
simulated  detonation  of  an  8-kT  source  as  computed  with  the 
SKIPPER  code. 


2.5  PARAMETRIC  STUDIES  OF  TUFF/WATER  MIXTURES 


A series  of  spherically  symmetric  one-dimensional 
shock  wave  propagation  problems  in  tuff/water  mixtures  were 
conducted  using  the  S3  SKIPPER  Lagrangian  finite  difference 
code.  These  calculations  were  made  to  evaluate  the  effects 
of  partial  and  full  saturation  of  the  tuff  matrix  on  the 
ground  motion  associated  with  an  underground  detonation  of 
8 kT  yield.  They  also  demonstrate  the  utility  of  the 
TAMEOS  routine  in  studies  of  this  kind. 


For  this  series,  the  spherical  SKIPPER  grid  was 
divided  into  400  zones  with  a maximum  radius  of  497.7  m. 
Thickness  of  the  first  zone  was  50  cm.  The  calculations 
were  carried  out  to  times  of  24  msec.  In  this  time  period, 
the  peak  stress  varied  from  an  initial  value  of  621  kbar 
down  to  less  than  1 kbar. 


Two  choices  for  total  pore  volume  f raction  . e . , ^n^  + ^j 
were  evaluated,  one  with  a total  porosity  of  20  percent  and 
the  other  with  a total  porosity  of  5 percent.  Water  content 
was  varied  so  that  the  air-filled  voids  accounted  for  100%, 

50-6,  and  0-6  of  the  total  pore  volume.  The  mixture  equations 
of  state  were  derived  from  the  tuff  and  water  equations  of 
state  described  in  Section  2.4.1.  TAMEOS  was  utilized  to 
compute,  store,  and  retrieve  the  mixture  states. 


The  plastic  crush  pressure  was  set  at  pc  = 20  kbar. 
Elastic  crushing  was  not  considered  (a  = a , p =0).  The 
deviatoric  response  of  the  mixtures  was  modeled^y  the  simple 
von  Mises  condition  with  a constant  shear  yield  strength  of 
Yq  = 1.0  kbar  and  a rigidity  modulus  set  at  y = 45  kbar: 


S2  + S2  + S2  < 1 Y2 

1 2 3 “ 3 Q 


(2.35) 


where  ai  = -p  + Si  expresses  the  total  principal  stress  in 
terms  of  the  pressure  and  deviatoric  stress. 
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2.5.1  Source 


The  energy  source  in  these  tuff/water  calculations 
consists  of  a spherical  cavity  which  contains  a gas  obeying 
the  ideal  gas  equation  of  state 


(2.36) 


A value  of  y = 1.4,  corresponding  to  diatomic  gas,  was 
utilized. 

Initially,  the  radius  of  the  cavity  is  3.72  meters 
and  the  internal  energy  of  the  gas  is  33.5  * 10 1 * ergs/g,  an 
amount  of  energy  equivalent  to  the  yield  of  8 kT  of  high 
explosive.  The  flow  within  the  cavity  is  not  calculated. 
Rather,  it  is  assumed  that  the  pressure  within  the  gas,  and 
thus  the  stress  acting  on  the  cavity  wall,  is  uniform  during 
each  time  step  in  the  calculation.  Moreover,  the  gas  is 
assumed  to  undergo  an  isentropic  expansion  (or  compression) 
i . e . , 


(2.37) 


where  p(t)  and  R(t)  are  the  cavity  pressure  and  radius, 

respectively.  The  initial  value  of  3.72  m for  R was 

selected  because  it  approximately  represents  the  volume  of 

rock  vaporized  by  an  energy  release  of  8-kT  of  explosive. 

Initial  pressure,  p , was  621  kbar. 

0 

2.5.2  Peak  Radial  Stress 

The  variation  of  peak  radial  stress  with  distance  from 
the  cavity  center  is  shown  plotted  for  the  porous  tuff 
calculations  in  Fig.  2.8.  The  201  porosity  cases  are  given 
in  Fig.  2.9. 
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Fig.  2.9~-Peak  radial  stress  decay  with  distance 
from  source  (8  kT)  for  205  porous  tuff,  three 
degrees  of  water  saturation.  Also  shown  are 
points  for  fully  saturated  51  porous  tuff  (from 
Fig.  2.8). 
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At  close-in  radii,  the  partially  and  fully  saturated 
matrices  undergo  slightly  higher  stress  levels  than  the  dry 
matrices.  As  the  initial  wave  propagates  out  to  larger  radii, 
and  the  stress  falls  below  pc  = 20  kbar,  the  partially  saturated 
matrix  peak  stress  levels  are  closer  to  those  of  the  dry  porous 
matrix  than  to  those  in  the  saturated  material. 

Most  significant,  peak  stresses  in  the  fully  saturated 
mixtures  are  very  close  to  one  another  despite  the  big  difference 
in  porosity.  Note  that  values  from  the  peak  stress  curve  for 
the  fully  saturated,  SI  porous  tuff  have  been  plotted  in 
Fig.  2.9,  and  are  very  close  to  the  201,  fully  saturated  case. 

It  is  clear  from  these  results  that  porosity  effects 
are  most  pronounced  when  the  degree  of  water  saturation  is 
reduced.  The  radial  stress  wave  is  more  retarded  by  the 
energy  absorption  associated  with  the  air-filled  voids  than 
with  the  presence  of  water.  In  the  201  porous  tuff,  this 
results  in  about  a 25 -kbar  reduction  in  the  peak  stress  for 
the  dry  matrix  within  8 meters  of  the  cavity  center.  At 
larger  distances,  the  partially  saturated  and  dry  matrices 
exhibit  stress  levels  relatively  far  below  the  fully  saturated 
material.  (Roughly  speaking,  the  peak  stress  is  reduced  by 
about  70  percent  between  30  and  10  kbar.) 

One  may  conclude,  therefore,  that  the  volume  fraction 
°f  air-filled  pores  is  the  most  influential  parameter  on  the 
peak  stress  levels.  The  half-saturated  media  responds  half- 
way between  that  of  the  saturated  and  dry  media  at  high  stress 
levels,  p > pc,  but  veers  closer  to  that  of  the  dry  material 
at  the  lower  stress  levels. 

2.5.3  Cavity  Growth 

The  cavity  radius  as  a function  of  time  in  each  calcu- 
lation is  presented  in  Figs.  2.10  and  2.11.  In  this  time  period, 
the  cavities  grew  monotonically . For  each  of  the  two  tuff 
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porous  tuff.  Three  degrees  of  water  saturation  of  51  porous  tuff  are  considered. 
The  cavity  radius  is  originally  3.72  m. 


porosity  values,  cavity  growth  was  maximum  for  the  completely 
dry  matrix  and  minimum  for  the  saturated  tuff.  Approximately 
5 percent  bigger  cavity  radius  was  produced  in  the  5%  porous 
tuff  (Fig.  2.10),  while  a 103  larger  cavity  (radii)  were  ob- 
served for  the  203  porous  tuff  matrices  (Fig.  2.11).  ' 

These  results  are  to  be  expected  since  the  higher 
percentage  of  air-filled  voids  leads  to  greater  crushing  of 
the  matrix  at  the  high  stress  levels.  One  could  anticipate 
a non-linear  relation  between  relative  cavity  wall  displacement 
and  porosity  due  to  spherical  geometry  effects. 

Once  again  it  may  be  observed  that  the  single  most 
important  parameter  is  ^n\  the  volume  fraction  of  air-filled 
voids.  The  cavity  radius  for  the  saturated  53  porous  tuff 
media  is  plotted  for  comparison  to  the  203  porous  result  in 
Fig.  2.11.  There  is  very  little  difference  between  these  two 
curves  and  it  is  insignificant  in  comparison  to  the  effect 
due  to  increases  in  the  void  content  or  the  media. 

2.5.4  Stress  vs  Time  Profile  at  R = 40  m 


It  is  at  radial  distances  much  greater  than  t lie  cavity 
radius  that  the  largest  relative  d :rences  occur  in  the 
stress- time  histories  of  the  various  matrices.  Fig.  2.12  is 
a plot  of  radial  stress  as  a function  of  time  at  the  40-m 
station  for  the  53  porous  tuff  matrix.  The  stress  levels 
at  this  location  were  below  the  crush  pressure  of  p = 20  kbar, 
hence  the  crush  process  is  incomplete  in  all  cases. 

The  stress  wave  in  the  fully  satuiated  mixture  travels 
quicker  and  reaches  the  40  m location  0.7  msec  prior  to  the 
crush  wave  associated  with  the  two  matrices  with  voids. 
Apparently,  the  low  degree  of  porosity  leads  to  insignificant 
differences  in  shock  arrival  time  for  the  latter  two  examples. 
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Clearly,  the  effects  of  increased  void  content  are 
evident  in  the  reduced  peak  stress  and  slower  rise  times 
of  the  stress  wave.  uteres tingly , after  the  first  4 msec 
following  wave  arri  in  the  saturated  media,  it  is  diffi- 
cult to  discern  significant  differences  in  the  stress  levels 
of  the  three  different  matrices. 

In  Fig.  2.13  is  the  analogous  stress  time  history  plot 
for  the  higher  porosity  tuff  matrix  calculations.  The  effects 
noted  for  the  5%  porous  case  are  more  evident  for  a 20%  poro- 
sity. The  wave  is  more  effectively  attenuated  by  the  matrices 
with  voids,  especially  the  totally  dry  case.  Once  again,  be- 
yond a certain  time  (19  msec),  it  is  difficult  to  differentiate 
between  the  three  media  during  the  unloading  process. 

From  these  results,  it  appears  that  radial  stress  attenua 
tion  and  rise  times  to  the  peak  during  the  loading  process  are 
extremely  sensitive  to  the  volume  fraction  of  gas-filled  pores. 
The  shape  of  the  release  curve  from  the  peak  stress,  however, 
is  quite  insensitive.  It  should  also  be  noted  that  the  fully 
saturated  media  exhibit  similar  stress-time  histories  (see 
Fig.  2.13).  Air-filled  porosity  accounts  for  the  major 
differences  in  the  shape  of  the  stress-time  profiles. 

2.5.5  Ground  Motion  of  Media  at  R = 40  m 

The  attenuated  stress  waves  associated  with  less 
saturated  media  lead  to  much  smaller  radial  displacements. 

This  is  portrayed  in  the  radial  displacement  curves  in 
Figs.  2.14  and  2.15.  Note  that  the  two  saturated  materials 
exhibit  almost  identical  displacement  histories  (see  Fig. 

2.15).  This  corroborates  the  main  conclusion  that  differences 
in  water  content  do  not  radically  affect  ground  motion 
characteristics  when  the  tuff  is  fully  saturated.  Air-filled 
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Fig.  2.13- -Radial  stress  time  histories  in  zone  originally  at  R = 40  m for  20% 
porous  tuff;  three  degrees  of  water  saturation.  (Also  shown  are  stress-time 
points  for  fully  saturated,  51  porous  tuff.) 


porosity,  however,  does  dramatically  alter  the  stress  wave 
propagation  characteristics  of  tuff  and  should  be  carefully 
modeled  in  ground  motion  calculations.  In  Appendix  B the 
effect  of  air-filled  porosity  is  studied  for  a high-explosive 
energy  source  buried  in  tuff.  This  effect  has  also  been 
recently  studied  by  Bjork^21^  and  Anderson.  ^2  .23  J 


50 


III.  IMPROVED  PREDICTIVL  METHODS  FOR  GRANITE 

3.1  INTRODUCTION 

A single  constitutive  theory  that  combines  the  higii 
pressure  equation  of  state  for  rock  with  a flow  law  which 
is  valid  at  low  pressures  and  high  strains  is  described  in 
this  section.  It  nas  led  to  a computational  scheme  suitable 
for  calculating  spherical  explosion  phenomena  in  both  the 
near  and  far  'ield.  The  theoretical  background  for  this  de- 
velopment, described  in  detail  in  a topical  report  by 
Dienes  , is  summarized  in  Sections  3.2  and  3.4.  Section 

3.2  presents  a theory  of  finite  strain  suitable  for  describ- 
ing large  distortions  of  porous  materials.  Section  3.3  lists 
tne  equations  of  motion.  Section  3.4  includes  several  al- 
ternative models  for  determining  the  flow  stress  in  rock 
masses,  and  a description  of  the  higii  pressure  equation  of 
state  used.  Also  included  is  an  approacn  to  estimating  strain 
rate  and  temperature  effects.  The  important  considerations  made 
in  incorporating  the  theory  into  the  SKIPPER  code  are  described 
in  Section  3.5.  Available  ground  motion  measurements  for 
underground  shots  in  granite  are  reviewed  in  Section  3.6. 

The  results  of  several  calculations  are  presented  and  then 
compared  with  the  underground  shot  data  in  Section  3.7.  By 
adjustment  of  the  flow  stress  it  was  found  possible  to  cal- 
culate a cavity  of  the  correct  size  and,  for  one  of  the  flow 
models  studied,  the  resulting  displacement  histories  are  in 
good  agreement  with  measurements.  Specifically,  the  kinematic 
hardening  rule,  wnich  models  anisotropic  hardening,  leads  to 
a realistic  calculation  of  shot  data.  The  "cap  model,  which 
assumes  isotropic  hardening,  results  in  a computed  behavior 
which  does  not  exhibit  the  observed  overshoot  in  displace- 
ment, and  shows  too  rapid  an  attenuation  of  the  shock. 

Kinematic  hardening  results  in  lower  hysteresis  and  slower 
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relief  waves  than  isotropic  hardening,  and  these  effects  ap- 
pear to  be  sufficient  to  bring  the  calculated  results  into 
better  agreement  with  measurements.  In  Section  3.8  the 
results  of  seven  calculations  are  presented  to  illustrate 
the  sensitivity  of  ground  motion  to  granite  material 
parameters  and  plasticity  models.  Finally,  Section  3.9 
discusses  the  state  of  the  research  to  date  and  its  rela- 
tion to  other  approacnes. 
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3.2  FINITE  DEFORMATION  THEORY 


Definitions  of  strain  can  be  obtained  in  a variety  of 
ways.  For  small  deformations,  they  are  generally  equivalent, 
but  there  may  be  substantial  differences  at  large  distortions, 
and  it  is  important  in  explosion  work  to  select  a definition 
which  will  lead  to  credible  solutions  for  arbitrarily  large 
motions.  Such  a definition  is  obtained  if  rate  of  strain  is 
defined  as  the  symmetric  part  of  the  velocity  gradient,  which 
in  tensor  notation  is  expressed  as 


+ u 


J 


(3.1) 


where  u^  denotes  the  ith  component  of  the  velocity  vector 
and  u.  . its  spatial  derivative  with  respect  to  the  coor- 
dinate  x^ . An  equivalent  definition  involves  representing 
the  velocity  gradient  as  the  sum  of  symmetric  and  antisymmetric 
parts , 
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(3.2) 


where  is  the  spin  tenror.  Its  components  have  the  same 

magnitude  as  those  of  the  vorticity  vector.  The  tensor,  D, 
whose  components  are  D^j  has  also  been  termed  the  rate  of 
deformation  tensor  and  the  stretching . ^ ^ 


In  spherical  symmetry,  the  case  of  current  interest, 
the  velocity  gradient  can  be  written 
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(3.3) 
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where  v is  the  radial  velocity.  To  obtain  explicit  ex- 
pressions for  the  strain,  we  may  adopt  the  Lagrangian 
description  of  the  deformation,  thereby  following  particle 
paths.  Then  the  previous  equation  reduces  to 
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where  r(r  ,t)  is  the  radius  of  a particle  initially  at  r . 

o 0 

If  the  definition  of  strain  for  spherically  symmetric  flow  is 

taken  to  be 
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it  can  be  readily  shown  that  the  relation 


e . . = D.  . 


(3.6) 


between  strain  rate  and  the  stretching  is  satisfied.  This 
equation  is  not  satisfied  in  alternative  treatments  of  finite 
deformation  plasticity.  For  example,  in  the  approach  followed 
by  Clifton1  J the  principal  strains  are  defined  by 


where  the  iL  are  the  principal  values  of  the  stretch  tensor. 
In  view  of  the  kinematic  identity 

D = j R (U  IT1  + U'1  U)RT  (3.7) 

given  by  Truesdell  in  Ref.  26,  the  linear  definition  of  strain 
given  above  is  not  consistent  with  Eq.  (3.7).  [In  flows  with- 
out rotation,  the  rotation  matrix  R reduces  to  the  identity 
matrix.  This  is,  of  course,  the  case  for  the  spherical  flows 
of  interest  here.] 


Although  in  spherical  flows  materials  may  undergo  large 
shears,  there  is  no  rotation  of  the  elements.  This  makes  it 
possible  to  separate  the  deformation  into  an  elastic  and  a 
plastic  part  in  a simple  way.  In  view  of  the  symmetry,  the 
stretch  tensor  has  the  simple  written  form 
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In  Ref.  24  it  is  shown  that 
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As  is  common  in  plasticity,  it  is  assumed  that  the 
strain  can  be  represented  as  the  sum  of  an  elastic  and  a 
plastic  part 
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(3.10) 


Expressing  the  elastic  and  plastic  parts  of  the  strain  as 
logarithms  of  the  corresponding  stretches, 


we  find 
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It  is  shown  in  Ref.  24  that  the  rate  of  change  of  compression, 
0,  can  be  expressed  in  the  form 

® ’ p/p  = ■(",  + F2  + F3)  • (3.14) 

Consequently,  the  compression  can  be  expressed  as  the  sum  of 

an  elastic  and  a plastic  part,  as  are  the  individual  strains, 
so  that 
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(3.15) 
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(3.16) 

The  plastic  part  of  the  deformation  of  a compact  material  does 
not  normally  involve  a volume  change.  Here  we  will  consider  the 
matrix  of  a porous  material  as  compact,  and  as  a result  the 
plastic  part  of  the  compression  can  be  interpreted  as  the  change 
in  void  volume.  The  void  fraction  is  defined  as  the  volume  of 
voids  per  urit  volume  of  material,  and  can  be  expressed  in 
terms  of  the  plastic  stretches  by 
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vv^ere  pu  initial  density  of  the  porous  material  and 

Pq  is  the  initial  density  of  the  matrix  portion  of  the  porous 
material . 
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3.3  EQUATIONS  OF  MOTION 


For  the  case  of  point  symmetry,  which  is  appropriate 
for  spherical  explosions , the  equations  of  motion  are  well 
known.  They  are  given,  for  example,  by  Wilkins or  in 
Ref. 24  . The  equation  of  mass  conservation  discussed  in 
the  previous  section  is  equivalent  to  the  rate  form 


p _ 2v  . 3v 
p r Tr 


(3.18) 


given  by  Wilkins  in  which  v is  the  radial  velocity.  The 
momentum  equation  is 


pv 


(3.19) 


and  the  energy  equation  is 


= °r  & + 2oe  7 


(3.20) 


which  is  equivalent  to  the  first  law  of  thermodynamics. 
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3 ' 4 CONSTITUTIVE  EQUATIONS  FOR  STRAIN  HARDENING  MATE RIALS 


In  spite  of  much  experience  witn  geologic  materials, 
the  formulation  of  a general  theory  suitable  for  deter- 
mining the  response  of  rock  or  soil  to  an  arbitrary  load 
has  received  intensive  study  only  recently.  The  problem  is 
more  complex  than  for  metals,  in  which  plastic  flow  takes 
place  at  nearly  constant  shear  stress.  In  rocks  and  soils  the 
shoar  strength  varies  dramatically  with  mean  stress,  and  the 
presence  of  voids,  cracks,  faults  and  pore  water  further  com- 
plicates the  material  description.  The  approach  taken  in 
this  research  was  to  modify  plasticity  theory  to  account  for 
these  complications  as  they  are  required.  In  many  respects 
the  theory  is  conceptually  more  elaborate  than  is  customary 
in  rock  mechanics.  The  reason  for  this  is  that  finite  strain 
and  high  energy  effects  are  accounted  for,  as  well  as  compaction 
and  dilatancy.  These  effects  are  believed  to  be  of  great 
importance  in  calculating  the  consequences  of  underground 
nuclear  explosions. 

The  approach  to  deriving  a sufficiently  general  consti- 
tutive relation  to  account  for  all  these  effects  draws  on  the 
thermodynamic  equation  of  state  for  the  high  pressure  behavior 
and  on  rock  mechanics  for  the  response  at  low  pressure. 

Several  theories  for  the  hardening  behavior  of  rocks  were 
studied  as  alternative  approaches.  These  theories  have  in 
common  that  they  provide  an  equation  for  the  yield  surface 
which  depends  on  the  history  of  the  deformation.  In  all 
cases  the  flow  law  of  von  Mises  was  used  in  connection  with 
the  equation  for  the  yield  surface  to  arrive  at  a constitutive 
relation.  This  relation  is  completed  by  using  the  constraint 
on  the  stresses  that  they  must  lie  either  on  the  yield  surface 
or  in  its  interior. 
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3-4. 1 Equation  of  State 


In  the  condensed  region  it  is  assumed  that  the  equation 
of  state  has  a form  similar  to  that  given  bv  Allen for 
geologic  materials,  which  expresses  the  pressure  as 

F = F(E,p)  = GEp  + f(p)  (3.21) 

where  p designates  the  matrix  density,  E designates  the 
specific  internal  energy  and 

G = a + ~K""b",  * n - £-  . (3.22) 

E n2  ° 

o 

In  Allen's  model,  which  emphasizes  the  fit  to  high  pressure 
data 


f(p)  " A p + Bp 2 , p = E i (3.23] 

^ Q 

To  match  the  low  pressure  behavior,  where  porosity 
introduces  a softening  effect,  Sandler  and  DiHaggio [ 301  use 
a variable  bulk  modulus 

A ■ \ l1  - a,e  0 ‘)  • (3.24) 

This  relation  can  be  put  into  equivalent  form  specifying  the 
pressure  as  a function  of  compression  if  we  observe  that 

Ji  = ' 3p  and  A = 3p"  • Then  a straightforward  integration 
leads  to 
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An  equation  for 

t(p)  which  fits  both  the 

low  and 

high  pressure  behavior 

is  given  by 

, v 3B  A 9\ 

fCp)  . ln(. 

0 

+ ( 1-a  )e  0 0 1 + B82  . 

(3.26) 

Here,  8,  the  matrix  compression,  is  given  by  Eq.  (3.15),  the 
superscript  being  temporarily  dropped.  The  last  term,  which 
is  quadratic  in  0,  should  be  dropped  for  underdense  states 
since  it  would  lead  to  unrealistic  behavior.  In  this  equa- 
tion, the  logarithmic  compression,  9,  rather  than  the  linear 
compression,  u,  is  now  used  to  allow  a treatment  compatible 
with  finite  deformation  theory. 


3.4.2  Hardening  Models  for  Rock  Behavior 


In  the  theory  of  plasticity  it  is  assumed  that  a 
definite  surface  exists  in  stress  space  with  the  property  that 
the  stresses  must  lie  either  on  the  surface  itself  or  in  its 
interior.  For  states  of  stress  represented  by  points  inside 
the  yield  surface,  the  material  is  described  by  the  theory  of 
elasticity,  or  possibly  a generalization  which  involves  non- 
linear elastic  and  thermo-elastic  effects.  For  states  of  stress 
on  the  yield  surface,  the  flow  is  governed  by  a flow  law  which 
constrains  the  stress  to  lie  on  the  yield  surface  if  tne 
elastic  theory  would  take  the  stress  outside.  The  flow  law 
of  von  Mises 


(3.27) 


leads  to  a unique  solution  in  most  cases,  with  X 
undetermined  multiplier  which  is  specified  by  the 


being  an 
cons  traint 
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When  the  function  specifying  the  constraint  is  not  equivalent 
to  the  function  f appearing  in  the  flow  rule,  the  flow  rule 
is  said  to  be  non-associated,  but  only  the  case  where  f = g 
was  investigated  in  the  current  study. 

The  increase  in  shear  strength  at  large  strains  that 
is  observed  in  most  materials  can  be  modeled  by  letting  the 
strength  depend  on  the  plastic  work,  the  second  invariant  of 
plastic  strain,  or  possibly  one  of  the  other  parameters  charac- 
terizing the  distortion.  If  it  is  assumed  that  the  flow  does 
not  depend  on  the  third  stress  invariant,  the  yield  surface 
has  the  hydrostat  as  an  axis  of  symmetry.  Models  have  been 
proposed  in  which  the  third  invariant  plays  a role,  such  as 
the  one  discuss''’  by  Cherry  ^31^  for  rocks  and  Freudenthal  ^ 32^ 
for  metals,  but  in  this  report  the  emphasis  is  on  yield 
surfaces  symmetric  about  a straight  axis.  Three  cases  are 
treated  in  the  discussions  that  follow.  In  the  first,  iso- 
tropic work  hardening,  the  yield  surface  is  a cone  which 
expands  isotropically  as  a function  of  the  plastic  work  done. 

In  the  second,  the  "cap"  model,  a portion  of  the  yield  surface 
is  conical,  but  it  is  completed  by  an  elliptical  cap.  The 
conical  portion  is  fixed,  but  the  cap  is  free  to  move,  though 
only  in  such  a way  that  the  flow  stresses  increase  in  magnitude 
Finally,  a kinematic  hardening  model  is  described  in  which  the 
conical  yield  surface  is  free  to  translate  in  stress  space  but 
not  to  deform.  These  models  are  conceptually  illustrated  in 
Fig.  3.1. 

The  stress  tensor,  a,  is  given  by  its  component  o— , 
and  it  is  often  convenient  to  write  the  governing  equations 
in  terms  of  the  components  rather  than  the  stress  tensor  itself 


(3.29) 


Defining  the  deviatoric  stress  by 


Sij  = °ij 


- J S../3 
i *3 


where 


is  the  first 


Hooke's  law' 


J,  * °u 


stress  invariant,  and  the  deviator  strain 

e . . = £••  - Ei.  i . 5 • • / 3 , 
ij  13  kk  13 

for  shear  deformation  is  expressed  as 


(3.30) 


by 

(3.31) 


where  u is  the  shear  modulus  of  the  material. 

When  the  shear  stress  attains  a critical  value  the 
flow  becomes  inelastic  and  is  then  governed  by  a constitutive 

law  of  the  form 


e? j = a 544  + ba-^  + 


*3 


13 


caik  akj 


(3.33) 


where  a,  b,  and  c are  usually  taken  to  be  functions  of  the 
stress  invariants  only.  Higher  order  terms  are  not  necessary, 
since  they  can  be  expressed  in  terms  of  the  lower  order  terms, 
according  to  the-  Cayley-Hamilton  theorem.  The  flow  la*  is 
further  restricted  if  we  adopt  the  form 


c?.  - X . C3.34) 

13  3oi3 

In  the  theory  of  metal  plasticity  it  is  assumed  that  = 0 

a condition  that  we  shall  call  isochoric  plasticity,  but  in 
an  investigation  of  the  stability  of  soil  masses  by  Drucker 
and  Prager^33]  this  constraint  was  lifted.  In  their  analysis 
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the  yield  surface  was  given  by 


f = >P7  ' gfJi)  = 0 (3.35) 

where 


J 

2 


(3.36) 


is  the  second  invariant,  and  the  Mohr-Coulomb  flow  condition 


) = Y - aJ 

1 i 


(3.37) 


was  adopted,  with  V the  yield  stress  in  simple  shear. 
Substitution  of  the  expression  given  above  for  f leads 
the  constitutive  equation 


to 


(3.38) 


With  this  result,  the  rate  of  change  of  specific  volume 
is  given  by 


(3.39) 


In  the  case  of  a Mohr-Coulomb  material 


g ' = -a 

where  a is  a positive  constant,  and 


(3.40) 


It  can  be  shown  that  1 is  always  positive,  and,  consequently, 
the  flow  always  exhibits  dilatancy  in  this  case.  The • dilatancy 
is  generally  greater  than  observed  in  tests,  and  in  reality 
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flows  in  which  compaction  occurs  are  also  observed.  This 
apparently  led  Drucker^-  ^ to  suggest  that  the  yield  sur- 
face should  be  completed  by  a spherical  cap,  for  in  that 
case  g can  have  either  sign.  At  high  pressures,  the  sign 
of  g'  is  positive  and  compaction  is  indicated.  At  low 
pressure,  g'  is  negative,  implying  dilatant  behavior.  This 
is  the  basis  for  assuming  that  porous  material  behavior  can 
be  described  by  the  methods  of  plasticity  theory. 

Sandler  and  DiMaggio  examined  in  detail  the  behavior 
of  granite  observed  by  Swanson,  and  proposed  an  analytic 
model  which  is  consistent  with  measurements.  The  model  is 
recapitulated  here,  since  it  represents  one  of  the  more 
promising  approaches  to  modeling  granite,  and  it  can  be 
conveniently  generalized.  The  yield  surface  is  represented 
in  two  parts,  a fixed  "failure  surface"  given  by 

8 * Yc  + (\  - Y le^1  (3.4  2) 


and  a variable  cap  described  by  an  ellipse  tangent  to  the 
failure  surface.  These  are  illustrated  in  Fig.  3.2,  which 
shows  the  yield  surface  in  a "reduced  stress  space".  In  this 
space  the  horizontal  axis  represents  distance  along  the 
hydrostat  and  the  vertical  axis  represents  the  radius  of 
the  yield  surface.  The  ellipse  tangent  to  the  failure 
surface  at  = J^p  is  given  by 


where  Jc  specifies  the  center  of  the  ellipse,  R denotes 
the  ratio  of  major  to  minor  axes  and  /$  is  the  length  of 
the  semi-minor  axis.  The  condition  that  the  cap  and  failure 
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Sketch  of  yield  surface  with  a movable  cap 
in  a "reduced  space  space". 


surface  be  tangent  at  their  intersection  determine  J and 

Q:  c 

Jc  " J i F " R2  8f  gF  (3.44) 

Q a Sp (i  + R2  gp2  ) C 3 . 45) 

where 
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-a  B Y e 1 lF  , 
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(3.47) 


and  is  the  value  of  at  the  point  of  tangency.  It 

was  determined  by  Sandler  and  DiMaggio  that  a good  fit  to  the 
data  can  be  obtained  by  defining  a hardening  parameter,  k,  bv 

* = U'gc)V^  £3-48) 

and  letting  the  abscissa  of  the  contact  point  move  according 
to  the  law 


J 


iF 


-Wk  . 


(3.49) 


For  the  shape  parameter 


they  assumed  a fit  of  the  form 


8 Wk 

R = R e 2 
o 


(3.50) 
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The  numerical  values  r,f  c j-. 

Cedar  City  Tonalite  are  ^ 

value  of  $ _nna  , l0w‘  (The  published 

*•  - tey°°  Ur-  *“ 

f°r  6 < > o.5  since  th  ^ modlfled  the  model 

surface  when  / takes  onY^  ^ th°  failure 

negative  hardening.  More  splcm^ny  T "T'  ^ t0 

rv;  * ■■  *•“  ••  w *** 

limiting  the  mavin,,,™  i 8‘  ‘5‘4  th  result  of 

shown.  as  dlscussed  above  is 


NUMERICAL  values  of  cap  model  parameters  for  cedar 
cm  TONALITE  GIVEN  BY  SANDLER  AND  DIMAGGIO 

[Ref.  3Q  J 


Y 

0 

152  ksi 

a 

i 

0.953 

R 

0 

4.0 

IV 

450 

8 

0 

0.002  ksi'1 

8 

i 

°-0029  ksi'1 

8 

2 

0.05  ksi'1 

u 

3300  ksi 

A 

0 

7500  ksi 

a 

0 

o 

L_ 

68 


u 

r) 


B 


Fig.  3.4--Cap  and  failure  surface  for  the  modified  cap  model  at  high 
stress  levels. 


3-4-3  Flow  Law  for  Materials  with  Hardening 

To  put  the  constitutive  equation  in  a form  suitable 
for  numerical  solution,  it  is  necessary  to  solve  for  the 
stress  rates  explicitly  in  terms  of  the  currer t stress  and 
the  strain  rate,  for  these  are  the  quantities  available  in 
the  computational  scheme.  This  is  accomplished  by  writing 
the  elastic  strain  rate  as  the  sum  of  a deviatoric  and  an 
isotropic  part 


V3  • 


(3. 51) 


The  elastic  shear  strain  is  given  by  Hooke's  law  in  rate 
fo  rm 


• U • 

= s../2p 


1J  “11 

and  the  matrix  volume  change  by 


xe  _ 3h  *T  . 3h  i 
6 ' ST-  J,  * If  E 


(3.52) 


(3.53) 


where  the  function  6e  = hGJ^E)  represents  the  equation  of 
state  of  the  matrix  material  in  a form  in  which  the  compression 
is  given  as  a function  of  J ^ , and  the  specific  internal 
energy,  E.  The  total  strain  rate 


+ 


can  then  be  written 


(3.54) 
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where 


T _ 1 . 1 3h 

T " 57  + T ST"  (3.56) 

i 

and 

J 

$ = s'  + Eg  • (3.57) 

To  complete  the  solution  it  is  assumed  that  A can  be 
written  in  the  form 


A = A,  + AJ1  (3.58) 

with  Aj  and  A not  explicitly  dependent  on  the  stress  rates. 
Then  the  constitutive  equation  can  be  expressed 

°ij  ‘ 2 °kk  ■ Aij  (3.59) 

where 


= 2u|eij  ■ \ 4^  *(\  ^ 


and 


♦ - 4 i)s 


1 9h 
I W 


(3.60) 


Zij  = 2yj(T  + AQ)  6ij  • Tg  aij}  (3-61) 

The  solution  of  Eq.  (3.59)  for  the  strain  rates  is 
straightforward,  following  the  usual  method  for  linear  simul- 
taneous equations.  For  spherical  explosions,  the  solution 
is  particularly  simple,  since  a = a , and  can  be  written, 

. A - 2R 

X = ii ,, 
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(3.63) 


where 


R = Z A - Z A (3.64) 

22  11  11  22 

and 

D = 1 - Z -21  . (3.65) 

11  2 2 

To  complete  the  analysis,  it  is  necessary  to  obtain  expressions 
for  A i and  A.  This  is  done  by  combining  the  results  of 
the  preceeding  analysis  with  the  equation  for  the  yield  sur- 
face, as  shown  in  the  sections  that  follow. 

3.4.4  Calculation  of  the  Multiplier 

The  flow  law  given  by  Eq.  (3.34)  leaves  an  undetermined 
multiplier,  X.  In  the  derivation  of  the  flow  law  based  on 
maximizing  the  rate  of  plastic  work  given  by  Hill  and 

due  to  von  Mises,  this  A appears  as  the  Lagrangian  multi- 
plier that  is  inevitably  introduced  in  extremum  problems 
involving  a constraint.  The  determination  of  A requires 
that  the  constraint  be  satisfied,  and  involves  differentiating 
the  constraint  equation  and  comparing  witn  an  appropriate  scalar 
equation  obtained  from  the  flow  law.  An  explicit  expression 
for  A ha  ing  the  form  A + AJ  Wu^  obtained  in  each  of 

i i 

the  four  cases  examined.  In  the  first,  the  constraint  is 
satisfied  in  the  absence  of  hardening.  The  remaining  three 
cases  involve  hardening  models,  and  are  simple  work  hardening, 
the  "cap"  model  and  kinematic  hardening. 

3 . 4 . 4 . 1 Multiplier  for  Materials  Without  Hardening 

To  obtain  an  expression  for  the  multiplier,  A,  the 
relation  obtained  by  differentiating  the  yield  condition, 
f = 0,  with  respect  to  time, 

J2  = 2gg ' jj  (3.66) 
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is  combined  with  a reduced  form  of  the  flow  law  obtained  by 
multiplying  Eq.  (3.55)  by  a-jj>  resulting  in  the  expression 


n 4.  J1  ahU 

p * 7-  -^)E  - 

1 9h  \ 

J JT 

g 

- S'  J. 

(3.67) 


which  has  the  form  X ^ + AJ  previously  assumed. 

3.4.4. 2 Multiplier  for  Materials  with  Isotropic 
Hardening 

If  the  yield  surface  is  assumed  to  expand  isotropically 
as  a function  of  the  work,  Wp , done  against  the  plastic  stresses, 
it  is  possible  to  represent  the  hardening  of  materials  in  a 
simple  fashion.  It  is  appropriate  to  note  here  that  the  un- 
loading behavior  is  not  well  represented  for  many  materials 
with  this  model,  since  the  Bauschinger  effect  is  not  accounted 
for.  Another  limitation  is  that  for  materials  whose  strength 
depends  significantly  on  the  mean  stress,  this  model  exhibits 
too  much  dilatancy , but  it  has  the  advantage  of  being  simple 
to  deal  with.  The  plastic  work  is  determined  by 


pWP  = o • . e? . 

ij  il 

and  the  yield  surface  has  the  equation 


(3.68) 


V7!  = *(v  wP) 


(3.69) 


Differentiating  this  relation  leads  to 


JL-  A- j + wp 

1 awP 


(3  70) 
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and  from  Eq.  (3.38),  we  can  show  that 


pWP  = \(g  - §f-  jJ.  (3.71) 

1 

The  rate  of  plastic  work  can  be  eliminated  from  these  equa- 
tions, and  after  some  straightforward  algebraic  manipulations 
a form  for  X 


of  the  required  type  is  obtained. 


(3.72) 


3 . 4 . 4 . 3 Multiplier  for  the  Granite  Cap  Model 

The  yield  surface  for  granite  formulated  by  Sandler 
and  DiMaggio  involves  both  a fixed  portion  and  a variable 
cap.  The  latter  is  represented  by  a family  of  ellipses  that 
varies  with  the  hardening  parameter,  ic,  given  by  Eq.  (3.48). 
When  the  stress  lies  on  the  fixed  portion  of  the  yield  sur- 
face the  analysis  for  flow  without  hardening  discussed  in 
Section  3. 4. 4.1  applies.  On  the  cap,  the  strerses  satisfy 
the  equation 


J 

2 


cv-o 


(3.73) 


where  gc  is  given  by  Eq.  (3.43).  It  is  easily  shown  that 
the  plastic  flow  law  of  Eq.  (3.27)  leads  to 


(3.74) 


75 


k 


To  determine  X the  flow  equation  is  operated  on  with  , 

the  time  derivative  of  Eq.  (3.73)  is  taken  and  < is  eliminated 
from  the  resulting  equations,  leading  to  the  result, 


X = 


i J J i 9h 

8C  3fic\*T 

(^p  T~  7£jl 

M 9J  ) 
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wc  gc  sgc  r pg 77  ’ 

y i + 3(^r) 


(3.75) 


3. 4.4.4  Multiplier  for  Kinematic  Work  Hardening 

r 3 7 1 

Prager's  original  rule1  J for  the  description  of  in- 
elastic behavior  assumed  that  the  yield  surface  translates 
in  stress  space  without  change  of  shape  as  plastic  deformation 
proceeds.  This  generalizes  the  bilinear  hysteresis  model 
sometimes  assumed  in  representing  plastic  deformation  when 
the  state  of  stress  is  one-dimensional.  In  the  current  ver- 
sion of  the  model  it  is  assumed  that  the  translation  of  the 
yield  surface  is  normal  to  its  axis  of  symmetry.  It  is  con- 
ceptually useful  to  visualize  the  surface  as  given  by  a 
vector  in  the  three-dimensional  space  of  principal  stresses. 
In  doing  this,  we  may  retain  the  double  subscript  notation, 
though  in  principal  axes  only  the  diagonal  terms  are  non-zero. 
A point  on  the  axis  of  symmetry  of  the  yield  surface  is  re- 
presented by 


a 


ij 


(3.76) 


where  the  term  a — 
and  the  second  term 
The  total  stress  is 


represents  the  displacement  of  the  axis 
represents  the  distance  along  the  axis, 
given  by 


aij 


(3.77) 
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with  representing  the  distance  in  stress  space  from 

the  axis  of  symmetry  of  the  yield  surface  to  the  stress 
point.  The  yield  surface  is  specified  by 

f - g f J x ) = 0 (3.78) 


where 


(3.79) 


The  crucial  physical  assumption  is  that  is  linearly 

related  to  the  plastic  strain  rate  by 


(3.80) 


where  b is  a hardening  parameter  and  e?.  is  the  deviatoric 

• D ^ J — 

plastic  strain  rate.  Since  e*jj  is  a deviator,  a^=  0,  im- 
plying that  the  translation  of  the  yield  surface  takes  place 
normal  to  its  axis  of  symmetry.  To  evaluate  X the  flow  law, 
Eq.  (3.55),  is  operated  on  with  , and  after  straight- 
forward manipulations  we  are  lead  to  the  result 


(3.81) 

3.4.5  Temperature  and  Strain  Rate  Effects 

The  magnitude  of  the  flow  stress  depends  on  both  the 
state  of  stress  and  the  history  of  the  motion,  as  discussed 
in  the  preceding  section.  In  addition  it  depends  on  tempera- 
ture and  strain  rate.  The  justification  for  treating 
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temperature  and  strain  effects  together  lies  in  the  physical 
basis  for  the  strain  rate  effect,  which  is  simply  that  when 

t“e  energ y of  a SrouP  of  at°™s  exceeds  the  activation  energy, 
some  slip  takes  place  through  a rearrangement  of  atoms.  The 
application  of  this  Arrhenius  activation  energy  concept  is 
discussed  by  Handin  and  Serdengecti  and  Boozer for 
rocks  and  applied  in  detail  to  metal  deformation  by  Zhurkov 
and  Sanfirova  and  more  recently  by  Samanta.[41J  A com- 
plete justification  for  the  factor  on  theoretical  ground* 
for  either  rocks  or  metals  has,  unfortunately,  not  been  es- 
tablished, but  Samanta  discusses  formulas  of  the  type 


e = ve  X2-kTA- 


(3.82) 


with  v the  "activation  volume"  and  AH  the  "activation 
enthalpy"  with  the  parameters  depending  on  the  flow.  If  we 
identify  the  "flow  stress,"  a,  with  the  yield  stress  of  the 
earlier  sections,  then  an  expression  for  the  yield  stress  in 

terms  of  temperature,  strain  rate  and  two  material  constants 
is  obtained 


U 1 *11  L,  /V 

2 


t J»  O JJ 


in  which  E'  is  the  second  invariant  of  the  deviatoric  strain 
rate  tensor.  The  substitution  of  yfF'  for  i generalizes  the 
one-dimensional  approximation  to  multi-dimensional  flows,  and 
t ie  constants,  a,  v and  Y have  to  be  obtained  by  experi- 
ments, theoretical  arguments,  or  estimates  based  on  experience 
with  other  materials.  Some  guidance  is  obtained  from  the  fact 
that  V is  of  the  order  of  magnitu  'e  of  the  atomic  vibration 
frequency  for  simple  materials,  such  as  aluminum,  and  a » k/MI 
is  determined  for  simple  materials  by  the  empirical  result  that 
the  activation  enthalpy  is  very  nearly  the  sublimation  energy. 
This  is  the  case  for  the  materials  investigated  by  Zhurkov  and 
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Sanfirova  and  by  Samunta.  The  situation  is  more  complex  for 
geologic  materials  tl  an  for  metals  and  alloys  but  an  under- 
standing of  some  of  the  simpler  examples  is  helpful.  Strain 
rate  effects  are  important  in  calculating  the  effects  of 
underground  nuclear  explosions  for  two  reasons.  First,  in 
the  near  field  the  temperature  is  high  because  of  shock 
heating  and  lowers  the  strength.  Second,  the  time  scale  of 
material  motion  is  generally  much  longer  for  nuclear  blasts 
than  for  laboratory  tests.  Consequently  it  seems  important 
to  estimate,  at  least  roughly,  the  magnitude  of  strain  rate 
effects  and,  in  th*.  absence  of  specific  data  for  the  case 
of  interest,  it  is  necessary  to  use  an  approximate  relation- 
ship. In  the  granite  calculations  a value  of  a of  0.5  x 10‘4 
was  used.  This  is  probably  somewhat  high,  but  it  was  selected 
to  minimize  the  effect  of  strength  near  the  cavity  at  a very 
early  time  in  the  development  of  our  techniques.  Experiments 
to  determine  a realistic  value  of  cc  would  be  desirable. 

To  use  the  correction  factor  1 + c.T  An  #T/v  u ls 
necessary  to  determine  the  temperature  in  the  course  of  the 
calculation.  To  simplify  the  calculation,  the  effect  of  the 
deviatoric  stresses  on  heating  is  ignored,  and  only  the 
thermodynamic  variables  pressure  and  internal  energy  are 
accounted  for.  This  is  a good  approximation  where  the  tempera- 
ture is  high  and  its  influence  is  significant.  Where  the  in- 
fluence of  the  deviatoric  stresses  is  important  the  temperature 
is  low,  and  it  is  not  essential  to  make  an  accurate  estimate 
of  its  deviation  from  normal. 

A differential  equation  for  temperature  can  be  obtained 
by  thermodynamic  considerations.  For  any  internal  energy 
E(V,T) 

dE  = CydT  + (|y)  <1V  . (3.84) 
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Combining  with  the  thermodynamic  identity 


(lv)T  = T(l?)v  ' p (3,85) 

it  follows  that 

CydT  = dE  - - p J dV  . (3.86) 

A convenient  expression  for  dT  can  be  obtained  by  using  the 
thermodynamic  identity 


and  the  definition 


of  Gruneisen's  ratio,  resulting  in  the  expression 
which  is  used  to  update  the  temperature. 


(3.87) 


(3.88) 


(3.89) 


For  the  form  of  the  equation  of  state  given  by  Eqs. 
(3.2])  and  (3.22)  we  find 


r 


a 


(3.90) 


The  differentials  required  in  Eq. 


(3.89)  are  given  in  terms 
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of  calculated  quantities  and  the  time  step,  dt , by 


dE  = [(oM-q)eM  + 2(oi2-q)e22]  dt/p  (3.91) 

dV  = (°n+  2a22)dt/p  (3.92) 


where  q represents  a viscous  stress  (artificial  viscosity). 
These  equations,  together  with  the  previous  equation  for  dT, 
are  continuously  used  in  the  computer  program  to  update 
temperature . 
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3*5  SET-UP  OF  SPHERICAL  EXPLOSION’  CALCULATIONS  IN  GRANITE 


A number  of  calculations  were  carried  out  to  determine 
how  well  the  theoretical  framework  described  in  the  previous 
sections  predicts  the  observed  motion  of  granite.  A yield  of 
1 kT  was  selected  for  these  runs,  and  in  comparing  with  nuclear 
shot  data  the  yield  is  scaled  down  to  1 kT  by  the  standard 
cube-root"  scaling  law'.  The  initial  cavity  radius  has  been 
varied,  but  in  the  calculations  cited  here,  it  was  held  at 
1.5  meters,  which  gives  a volume  that  is  realistic  for  medium 
yield  shots. 

The  source  was  represented  by  a polytropic  gas  law 
with  an  index  y = 4/3.  This  value  was  selected  as  being 
representative  since  it  is  exact  for  pure  radiation,  which 
governs  for  very  high  yields  at  early  times, and  also  for  a s 
tri-atomic  gas  without  internal  degrees  of  freedom.  Since 
the  cavity  contains  large  amounts  of  water  and  silicon  dioxide, 
it  was  felt  that  the  real  values  might  not  vary  significantly 
from  4/3.  Studies  by  Allen  and  I)ufft42]  and  Wagner  and  Louie  [43J 
have  indicated  that  the  results  are  not  highly  sensitive  to  the 
value  of  y,  but  a precise  calculation  should  allow  for  ioni- 
zation and  other  non- ideal  features  of  the  cavity  gases. 

A source  calculation  may  involve  either  the  expansion 
of  thin  spherical  shells  of  gas,  or  it  may  assume  that  the  shells 
mix  as  the  result  of  turbulence.  If  the  mixing  is  thorough, 
it  is  a good  approximation  to  represent  the  source  as  a uni- 
form sphere  expanding  adiabatically , an.  this  is  the  approxi- 
mation adopted  in  the  current  study.  In  view  of  the  difficulty 
of  studying  mixing  in  a transient  flow  theoretically,  a descrip- 
tion of  the  source  Dehavior  will  probably  require  an  experimental 
approach.  The  discussions  by  Butkovitch  ^2f)  ’ 44  ^ indicate  that 
the  water  content  of  the  rock  significantly  affects  the  results, 
but  the  emphasis  in  these  calculations  was  on  the  effect  of 
rock  strength,  and  water  content  was  not  accounted  for. 
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The  zoning  for  the  1-kT  problems  is  illustrated  in 
Fig.  3.5.  The  computer  program  was  formulated  to  allow  1000 
zones,  but  only  350  were  required  for  these  calculations. 

The  first  zone  is  50  cm  in  thickness,  and  the  thickness  was 
increased  by  1 percent  per  zone.  The  total  radius  accounted 
for  in  the  calculational  grid  is 


(3.93) 


where  a is  the  thickness  ratio,  r is  the  initial  radius 

o 

of  the  source  region,  Ar  is  the  width  of  the  first  zone, 

o 

and  N is  the  number  of  zones.  This  radius  is  1563  meters 
for  the  indicated  zoning. 

The  SKIPPER  program  uses  the  standard,  Lagrangian 

approach  to  the  calculation  of  spherical  motions  described 
F2  8 1 

by  Wilkins.1  1 In  addition  to  the  stresses  resisting  the 
motion  described  in  the  preceding  sections,  an  artificial 
viscosity  is  included  in  the  program  which  incorporates  both 
linear  and  quadratic  terms.  Details  are  given  in  the  RIP 
report  by  Fisher,  Cecil  and  Lane^4^  (the  SKIPPER  code  origi- 
nated from  the  RIP  codt ) . It  was  found  in  the  course  of  the 
current  study  that  an  improvement  in  the  behavior  at  the  shock 
frojiZ  can  be  obtained  by  eliminating  only  the  quadratic  term 
in  the  artificial  viscosity 

q = (CQ  Au2  - CLcAu)p  (3.94', 


where  Au  is  the  velocity  change  across  a zcne,  c is  thr 
sound  speed  and  p is  an  average  density,  when  the  material 
is  expanding.  In  the  usual  approaun,  both  terms  are  dropped 
when  material  is  expanding,  but  this  leads  to  an  excessively 
non-uniform  treatment.  A comparison  of  the  two  methods 
Is  illustrated  in  Fig.  3.6,  which  shows  that  a mild  instability 
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-Initial  conditions  and  zoning. 


in  the  calculation  that  occurs  for  relatively  large  values 
of  artificial  viscosity  (C^  = 3.2,  Cj  = 1.0)  disappears 
when  the  linear  term  is  retained  in  the  expansion  region, 
even  though  the  viscosity  was  lowered  at  the  same  time. 
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3.6  UNDERGROUND  SHOT  DATA 


Before  entering  into  a discussion  of  the  calculated 
results,  it  seems  appropriate  to  discuss  some  of  the  data 
taken  from  the  Hardhat,  Shoal  and  Piledriver  shots,  since 
these  data  are  the  basis  for  evaluating  the  validity  of  the 
calculations.  This  is  especially  true  since  there  is  a long 
history  of  calculations  of  underground  shots  with  highly 
variable  degrees  of  success.  In  analyzing  the  shot  data  it 
is  observed  that  the  cavity  volume  per  kiloton  of  yield  is 
remarkably  consistent  between  shots.  The  data  are  summarized 
in  Table  3.2.  The  consistency  of  the  volume  per  unit  yield 

TABLE  3.2 

CRATER  RADIUS  AND  CRATER  VOLUME  PER  TON 
FOR  THREE  GRANITE  SHOTS 


Shot 

Yield, 

Kilotons 

Cavity 
Radius , 
Meters 

Reference 

Cavity  Volume 
Per  Unit  Yield, 
Cubic  Meters 
Per  Ton 

Hardhat 

5.0 

19.2 

46 

5.80 

Shoal 

12.5 

25.6 

46 

5.69 

Piledriver 

61. 

44.5 

47 

6.03 

far  exceeds  the  credibility  of  the  individual  measurements. 
Though  both  yield  and  cavity  volume  vary  by  nearly  201, 
depending  on  the  author  and  type  of  measurement,  the  average 
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cavity  volume  per  unit  yield  is  5.7  cubic  meters  per  ton 
with  a variation  of  only  4 percent.  For  a 1-kT  shot,  the 
radius  would  be  11.1  meters. 

[47  i 

Perret  in  his  discussion  of  Piledriver1  J and  Werth 
and  Herbst  in  their  discussion  of  Hardhat ^ ^ have  pointed 
out  that  the  late-time  displacement  can  be  estimated  from 
an  elementary  analysis  based  on  the  argument  that  the  mass 
of  material  inside  the  sphere  with  radius  R before  the 
shot  must  equal  the  mass  inside  the  sphere  of  radius  R+6 
after  the  shot,  where  5 is  the  displacement , provided  the 
the  material  does  not  undergo  significant  compaction  or 
bulking.  It  is  reasonable  for  the  granite  shots  cited  here 
to  assume  that  the  initial  cavity  volume  is  negligible.  The 
mass  equation  then  takes  the  form 

t tt(R+6)  3 - i ttR^  = ^ ttR  3 (3.95) 

which  reduces,  for  5 <<  R,  to 

6 = R(3/3R2  (3.96) 

where  R^  is  the  cavity  radius.  For  the  Piledriver  shot 
this  leads  to  a value  of  1.04  ft  for  the  permanent  displace- 
ment at  a radius  of  1000  ft.  The  fit  by  Borg^^  to  26  data 
points  results  in  a displacement  of  0.86  feet  at  this  radius, 
a discrepancy  of  about  201,  consistent  with  the  uncertainty 
in  displacement,  cavity  radius  and  yield.  The  need  for  a 
reliable  yardstick  such  as  this  in  correlating  theory  and 
shot  data  is  vital,  since  individual  measurements  may  vary 
by  a factor  of  five  from  the  mean  value.  Many  of  the  measure- 
ments of  displacement  cited  by  Borg,  vary  by  a factor  of  ten 
from  one  another.  Since  these  authors  of  granite  shot  studies 
agree  on  the  validity  of  the  above  displacement  formula,  it 
seems  to  provide  a good  benchmark  for  testing  the  validity 
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of  calculations.  In  view  of  the  uncertainty  in  measurements, 
it  is  felt  that  caution  should  be  exercised  in  searching  for 
a correlation  of  individual  measurements  and  calculations, 
such  as  velocity  history  at  a fixed  station. 

Displacement  history  may  be  one  of  the  more  reliable 
indicators.  Since  it  is  an  integrated  quantity,  it  washes 
out  some  spurious  details,  and  the  final  value  can  be  esti- 
mated from  the  formula  given  above.  This  eliminates  the 
possible  effect  of  slow  drift  in  the  traces  due  to  long  term 
electrical  effects. 

The  dynamic  overshoot  at  displacement  is  very  large. 
The  peak  displacement  is  compared  with  the  final  value  in 
Table  3.3,  and  an  overshoot  of  3501  is  apparently  typical. 

TABLE  3.3 

DISPLACEMENT  OVERSHOOT  FOR  THREE  GRANITE  MEASUREMENTS 


Piledriver 
Station  1 
(658  ft) 

Piledriver 
Station  2 
(1543  ft) 

Hardhat 
(457  m) 

Peak  Displacement,  d 

r max 

65  in. 

15.5  in. 

3.8  cm 

Final  Displacement,  dw 

13  in. 

9 in. 

(5.2  in.)* 

1.19  cm 

Overshoot  (100  dmax/dj 

4001 

172% 

(300%)* 

320% 

It  is  of  considerable  interest  to  explain  this  overshoot 
quantitatively  by  theoretical  methods,  but  the  constitutive  equa- 
tions previously  used  were  not  able  to  model  the  overshoot  and 
late-time  displacement  in  a straightforward  fashion  It  will  be 
shown  that  the  kinematic  hardening  model  does  succeed  in  providing 
a reasonable  theoretical  description  of  the  overshoot. 


The  9- in.  figure  given  is  questionable,  and  the  alternative 
value,  (5.2  in.)  based  on  Borg's  fit,  is  also  tabulated  for 
comparison. 
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3 . 7 COMPARISON  OF  CALCULATIONS  AND  I- 1 li  L D DA  TA 


Three  calculations  are  described  in  this  section  which 
bear  on  the  problem  ol  how  best  to  model  granite.  Two  of  these 
involve  the  same  material  constants,  where  applicable,  and  were 
designed  to  compare  the  "cap"  and  the  kinematic  hardening  models. 
In  these  cases  the  strength  (4.0  kbar,  reference  value)  w>as 
clios  en  so  that  the  "cap"  model  calculation  would  result  in  a 
cavity  of  roughly  the  right  volume.  Details  of  these  calcula- 
tions are  supplied  in  the  three  sections  that  follow.  The 
reference  value  for  strength,  Y^,  is  the  strength  in  the 
formula 

Y.  ’ hi1  * “ T *»#>) 

and  represents  the  strength  at  T = 0 and  J = 

i 

3.7.1  Cap  Model  Calculation 

From  previous  studies1"50^  and  test  runs  it  was  known 
that  the  parameters  given  in  the  original  Sandler  and  DiMaggio 
fit  to  Cedar  City  Tonalite  wrould  lead  to  an  unrealistically 
small  value  of  cavity  radius.  With  the  premise  in  mind  that 
the  cavity  radius  should  come  out  near  to  11.1  meters  for  a 
one-kiloton  shot,  as  discussed  in  Section  3.6,  the  yield 
strength  was  adjusted  to  give  a realistic  cavity  volume.  A 
set  of  parameters  that  accomplishes  this  is  listed  in  Table  5.4. 

In  the  table  is  also  included  a value  for  the  dilata- 
tional  wave  speed, 

= V(K  + 4p/5)/p  > K = Ao(l-a  ) . (5.97) 

The  value  of  aQ  was  adjusted  so  that  this  wave  speed  would 
be  consistent  with  the  4.8  km/sec  value  for  Mardhat  granite 
at  zero  pressure  published  by  Werth  and  IlerbstJ48^ 


For  the  choice  of  parameters  listed,  the  calculated 
crater  radius  is  10.5  meters.  This  is  close  enough  to  the 
11.1-meter  figure  of  Section  3 .6,  which  was  based  on  shot 
data,  to  justify  a detailed  comparison  of  other  measurements 
with  calculations.  The  main  observation  resulting  from  the 
comparison  was  that  the  cap  model  results  in  a displacement 
history  which  exhibits  virtually  no  overshoot,  whereas  the 
test  data  indicates  values  of  overshoot  exceeding  300%.  The 
detailed  results  of  the  calculations  are  deferred  to  Section 
3.7.3  in  which  a comparison  with  both  the  data  and  the 
kinematic  hardening  model  is  discussed. 

3.7.2  Kinematic  Hardening 

The  kinematic  hardening  model  discussed  in  Section 
3.4. 4.4  was  used  as  the  basis  for  a one-kiloton  calculation 
with  the  same  mechanical  parameters  as  those  used  in  the  pre- 
ceding discussion,  except  that  the  cap  behavior  was  not 
accounted  for.  Hardening  due  to  cap  motion  was  replaced  by 
translating  the  axis  of  symmetry  of  the  yield  surface  according 
to  Eq . (3.76)  with  displacement  given  by 

*ij  = beij  (3.98) 

the  value  for  b was  100  kilobars.  A slope  of  100  kbar 
is  shown  in  Fig.  3.7;  which  is  taken  from  Ref. 35  to  illustrate 
granite  behavior.  The  cavity  radius  and  displacement  showed 
substantial  overshoot  in  the  calculation,  and  the  final  cavity 
radius  was  smaller  than  w’ith  the  cap  model  of  hardening. 

It  was  estimated  that  a reduction  in  the  reference  strength 
from  4 to  1.68  kbar  would  increase  the  cavity  radius  to 
11.1  m.  The  result  of  these  calculations  are  included  in 
Figs.  3.8  through  3.13. 
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A.  Cycled  to  1/4  maximum  load 

B.  Cycled  to  1/2  maximum  load 

C.  Cycled  to  3/4  maximum  load 

D.  Cycled  to  1/4  maximum  load 


Fig.  3. 7- - 1 llustration  of  a stress-strain  relation 
taken  from  Ref. 35,  showing  the  Bauschinger  effect 
in  Cedar  City  tonalite,  and  the  theoretical  model 
for  pure  shear. 
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Fig.  3.8  Comparison  of  three  calculations  with  observed  cavity  size. 


-Displacements  at  51.9  m at  Station  1 of  Piledriver 
scaled  to  1-kT  yield,  based  on  data  of  Ref  47. 
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Displacement  at  120  m for  a yield  of  1 kT  and  comparison 
with  measurement  at  Station  2 in  Piledriver,  based  jn 
data  of  Ref.  47 


Fig.  3. 11- -Displacement  at  267  m from  a 1-kT  source  (scales  to  457 
from  Hardhat  source- - locat ion  of  measurement  3VR) , based  on 
data  of  Ref.  48 


Permanent  Radial  Displacement  (cm) 


Peak  Velocity  (ft/sec) 


Fig.  3. 13- - Comparison 
peak  velocities  for  a 


of  observed  and  calculated 
1-kT  explosion  in  granite. 


5.7.3  Comparisons  of  Calculations 


The  cap  mode]  can  be  adjusted  to  give  a realistic 

value  of  crater  size  by  lowering  the  value  of  flow  stress  from 

the  laboratory  value  of  10.56  kbar  to  2.6  kbar  at  normal 

strain  rates.  The  maximum,  or  reference,  flow  stress  is  4.0 

kbar,  but  the  2.6  kbar  value  is  cited  for  comparison,  since 

this  is  the  value  of  Y that  would  be  obtained  at  normal 

o 

conditions  using  the  strain-rate  relation  of  Lq.  (3.83).  (A 
justification  for  this  _s  that  large  masses  of  rock  can  be  ex- 
pected ' include  large  cracks  and  hence  a lower  mean  strength, 
than  laboratory  samples.  In  addition,  the  presence  of  water 
in  situ  can  be  expected  to  lower  the  strength  significantly.) 
The  time  history  of  crater  radius  using  the  cap  model,  with 
adjusted  strength  is  shown  in  Tig.  3.8.  The  final  value  of 
crater  radius  is  10.5  m,  whereas  the  best  estimate  based  on 
averaging  Piledriver,  Hardhat  and  Shoal  data  is  11.1  m.  In 
view  of  the  many  uncertainties,  it  did  not  seem  appropriate 
to  improve  this  agreement  by  adjusting  parameters.  In  the 
same  figure  the  history  of  cavity  radius  is  shown  for  kinematic 
hardening.  The  strength  is  unrhanged  from  the  cap  model  value, 
but  the  radius  is  somewhat  lowered,  and  a significant  amount 
of  overshoot  takes  place.  To  get  a computed  cavity  radius  that 
agrees  with  the  shot  data,  the  reference  strength  for  the 
kinematic  hardening  calculation  was  lowered  to  1.6S  kbar.  It 
can  be  seen  in  the  figure  that  the  amount  of  overshoot  was 
significantly  increased  as  a result  of  this  change.  Unfor- 
tunately, as  a result  of  the  way  the  problem  was  set  up  (the 
number  of  zones  was  limited  to  350) , it  was  not  possible  to 
continue  the  calculation  until  the  oscillation  died  out.  It 
seems  probable,  however,  that  the  final  value  of  cavity  radius 
would  be  near  the  observed  size. 

Comparisons  with  measured  displacements  are  shown  in 
Figs.  3.9,  3.10,  and  3.11.  The  first  two  correspond  to  the 
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668-  and  1543-foot  stations  for  Piledriver.  Th.  displace- 
ment histories  given  by  Perret  were  scaled  to  1 kT,  assuming 
a yield  of  61  kT  and  the  cube  root  scaling  law  to  get 
the  traces  caUed  "observed".  The  third  comparison  with  a 
measured  displacement  history  is  based  on  the  displacement 
trace  reported  by  Werth  and  Herbst.  In  the  first  comparison, 
at  a scaled  radius  of  51.9  m , the  "cap"  model  curve  approaches 
the  value  based  on  the  Borg  fit,  16  cm,  with  very  little  over- 
shoot, the  peak  amplitude  being  18  cm.  The  observed  peak  dis- 
placement, scaled,  is  41.5  cm.  Inspection  of  the  figure  shows 
that  the  kinematic  hardening  model  results  in  a more  realistic 
behavior.  The  peak  overshoot  is  50  cm,  for  = 1.68  kbar , 
and  the  timing  is  generally  like  that  observed.  At  a radius 
of  120  m the  conclusions  are  similar.  The  cap  model  results 
in  a peak  displacement  of  1.90  cm  where  the  measured  peak  is 
10.1  cm,  and  it  occurs  considerably  later.  The  kinematic 
hardening  model  gives  a peak  of  8.3  cm  and  the  timing  is  also 
roughly  similar  to  that  observed.  In  interpreting  these  re- 
sults, it  should  be  borne  in  mind  that  there  is  considerable 
variation  in  measurement  between  different  instruments  at  the 
same  station.  The  difference  between  the  kinematic  hardening 
calculation  and  the  observed  displacement  is  not  greater  than 
the  difference  between  records  from  a velocity  pick-up  and 
an  accelerometer  at  this  station. 

The  displacement  reported  by  Werth  and  Herbst,  scaled 
to  1 kT,  has  a peak  value  of  2.18  cm  where  the  cap  model  gives 
a peak  value  of  0.82  cm,  and  the  peak  occurs  somewhat  early. 
The  kinematic  hardening  model  gives  a peak  value  of  1.8  cm 
which  occurs  only  slightly  earlier  than  the  observed  peak 
value . 

The  permanent  displacement  data  is  summarized  in 
Fig.  3.12,  with  the  experimental  points  being  reduced  to 
1 kT  by  cube  root  scaling,  as  before.  The  values  for  Hardhat 
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and  Piledriver  displacement  are  taken  from  Perret's  report, 
except  for  the  last  Hardhat  point,  which  is  taken  from  Werth 
and  Herbst.  The  experimental  data  points  generally  correlate 
wel]  with  Borg's  fit  and  with  the  theoretical  value,  R3/3R2, 
except  for  one  Piledriver  point  which  is  questionable/  The' 
"cap"  model  predicts  more  attenuation  than  is  observed, 
whereas  the  kinematic  hardening  model  gives  very  nearly  the 
correct  slope.  Orly  the  kinematic  hardening  calculation  for 
the  reference  strength  of  4.0  kbar  is  plotted,  since  the  dis- 
placement had  not  reached  a steady  state  in  the  second  kine- 
matic hardening  calculation. 

Peak  velocity  data  are  compared  in  Fig.  3.13,  which  is 
taken  from  Perret's  report  and  compares  data  from  Piledriver 
Hardhat  and  Shoal,  scaled  to  1 kiloton.  Superimposed  on  the’ 
data  are  the  results  of  a cap  model  and  a kinematic  hardening 
model  calculation,  both  for  the  4.0  kbar  reference  strength. 

The  cap  model  calculation  lies  consistently  below  the  measure- 
ments, whereas  the  kinematic  hardening  model  is  in  general 
agreement  with  the  measurements.  The  slope  of  the  kinematic 
hardening  calculation  seems  to  be  a little  lower  than  the 
indicated  slope  of  the  measured  peak  velocities.  Unfortunately, 
peak  values  of  velocity  for  the  kinematic  hardening  model  with 
reference  strength  of  1.68  kbar  are  not  available  because  the 

printing  frequency  in  setting  up  the  calculation  was  too 
low. 


102 


3.8  PARAMETRIC  STUDIES  OF  GRANITE 


A series  of  spherically  symmetric  one-dimensional 
shock  wave  propagation  problems  in  granite  was  run  using 
the  SKIPPER  finite  difference  code.  The  purpose  of  these 
calculations  was  to  evaluate  the  effects  on  ground  motion 
of  varying  certain  material  properties.  To  run  problems 
in  which  all  parameters  occurring  in  the  constitutive  rela- 
tions are  varied  in  a controlled  manner  so  that  the  effect  of 
each  on  wave  propagation  could  be  evaluated  would  be  a very 
lengthy  task.  Inasmuch  as  the  greatest  uncertainty  in 
material  response  models  for  media  affected  by  an  explosion 
in  the  ground  is  for  stresses  in  the  region  below  ten  kilo- 
bars,  variations  of  strength  parameters  were  investigated. 
Thermodynamic  parameters,  which  would  be  more  important  in 
the  very  Wigh  pressure  regime,  were  held  constant.  Ideally, 
it  would  have  been  desirable  to  carry  the  problems  to  the 
point  at  which  the  peak  stress  levels  had  attenuated  to  the 
level  of  seismic  signals  and  motion  had  essentially  ceased. 
This  was  done  in  the  calculations  for  comparison  with  field 
data  (Section  3.7)  , but  was  not  feasible  for  the  parameter 
studies.  The  problems  were  run  to  times  of  about  20  msec 
and  stresses  had  attenuated  to  the  order  of  one  kilobar. 

These  calculations  were  sufficient  to  show  the  effect  of 
varying  certain  strength  parameters  on  peak  stress  attenuation 
cavity  growth,  pulse  shape  and  particle  displacement.  The 
strength  parameters  were  chosen  in  a range  that  is  represen- 
tative of  what  one  might  expect  to  find  in  granite. 

The  source  in  all  but  one  of  the  calculations  reported 
in  this  section  is  the  same  as  used  for  the  material  para- 
meter study  for  tuff  reported  in  Section  II.  It  consists  of 
a spheiical  cavity  in  the  material  which  contains  a gas 
obeying  a pressure-volume -energy  (p,V,E)  equation  of  state 


given  by 

E = (3.99) 

The  value  y = 1.4  was  used.  Initially  the  radius  of  the 
cavity  is  3.72  meters  and  the  internal  energy  of  the  gas  is 
33.5  x 10  19  ergs,  an  amount  of  energy  equivalent  to  the 
yield  of  8 kT  of  High  explosive.  The  flow  within  the  cavity 
gas  is  not  calculated.  The  pressure  within  the  gas,  and  thus 
the  stress  acting  on  the  cavity  wall,  is  taken  to  be  uniform 
during  each  time  step  in  the  calculation  and  is  given  by 

/ R n \3^ 

p(t)  - PjAirrh")  * p0  * 621  ,<bai  (3.ioo) 

where  p(t)  and  R(t)  are  the  cavity  pressure  and  radius 
rerpecitvely  at  time  t.  The  initial  value  of  5.72  m for 
the  cavity  was  chosen  on  the  basis  that  it  approximately  re- 
presents the  volume  of  rock  vaporized  by  an  energy  release  of 
8 kT  of  explosive.  In  one  problem  the  cavity  radius  was  in- 
creased to  8 m to  give  an  order  of  magnitude  increase  in  the 
volume  in  which  the  8 kT  was  deposited. 

The  p-V-E  equation  of  state  used  for  the  rock  material 
was  the  blend  between  the  high  pressure  and  low  pressure  forms 
described  in  Section  3.4.1.  The  required  material  constants 
used  are  those  listed  below: 

P = 2.68  g/cc 
o 

a = 0.5 
b = 1.3 

E = 1.6  x 1011  ervs/g 

o ■ ° 

A = 518  kbar 
o: 
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B = 180  kbar 

a = 0.409 
o 

3 = 0.029  kbar*1 

o 

Em  = 1,8  * 1Ql°  er8s/fi 

Except  for  a small  change  in  a and  A these  are  the 

0.  0 

same  values  used  for  granite  in  Section  3.7  (Table  3.4). 

A principal  stress  component  is  given  by 

oi  = -p  + Si  (3.101) 

where  is  the  associated  principal  deviatoric  stress. 

In  the  calculations  reported  in  this  section,  three  models 
were  used  to  limit  the  deviatoric  terms  and  treat  the  strength 
of  the  rock  material.  These  three  models  are: 

(a)  Simple  von  Mises 


S2 


S2 

2 


s2 

3 


if 


Y2 

o 


I'-U 


where  Y is  a constant  yield  strength. 
0 


(3.102) 


(b)  Mohr- Coulomb 


S2  + S2  + S2 

1 2 3 


(3.103) 


where  Y(p)  is  a yield  strength  which  is  dependent  upon  the 
hydrostatic  pressure.  In  the  present  work  the  form  of  Y(p) 
is  taken  to  be  the  exponential  form 


Y(p)  = Y - 
0 


where  Y , Y , and 
o 1 


8 are  constants. 

i 


(3.104) 
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1 


for  all  cases  the  rigidity  modulus  was  held  constant  at 
U = 228  Khars . The  factor  (1  - multiplying  the  Y's  rep- 

l esents  a thermal  softening  term.  When  the  internal  energy  L 
reaches  the  melt  energy,  F^,  the  strength  vanishes  and  the 
material  can  no  longer  support  any  deviatoric  stresses. 

(c)  Capped  Yield  Surface 

The  yield  surface  has  the  exponential  form  given  by 

Tq.  (3.104).  However,  it  is  capped  by  an  ellipse  which  is 

tangent  to  the  exponential  yield  surface  at  their  point  of 

intersection.  The  elliptical  cap  moves  in  the  Y-p  plane* 

according  to  the  amount  of  plastic  work  hardening  that  occurs. 

The  point  of  intersection  of  the  yield  surface  and  the  cap, 

the  ratio  of  the  major  and  minor  axes  of  the  elliptical  cap, 

and  the  center  of  the  ellipse  are  functions  of  a work  hardening 

parameter.  The  computation  of  the  work  hardening  parameter, 

k,  is  described  in  Section  3.4.2.  The  yield  surface  and  cap 

are  shown  in  Fig.  3.4. 



Some  confusion  can  arise  between  the  p-Y  plane  and  the 
Ji  vs  yJT  Plane  where  J1  is  the  first  stress  total  in- 
variant and  J2  is  the  second  deviatoric  stress  invariant. 
Here  the  following  definitions  are  used: 

P = - i ( o + o + a ) 

1 2 3 

4 y 2 = S2  + S2  + S2 

J 12  3 

J = a + o + a 

112  3 

J = i(s2  + s2  + s2) 

2 ^ ' 1 2 3 ' 

Thus,  J = -3p  and  \lj~  = — 

1 V 2 yfi 
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After  several  short  test  runs  to  assure  that  the  con- 
stitutive  relations  described  above  were  functioning  satis- 
factorily in  the  SKIPPER  code,  a total  of  seven  production  runs 
calculations  on  granite  were  carried  out  to  times  in  excess  of 
20  milliseconds.  A description  of  the  plasticity  models  for 
these  runs  is  summarized  in  Table  3.5  and  the  calculational 
results  are  presented  in  Figs.  3.14  through  3.25. 

The  asymptotic  value  approached  by  the  shear  yield 
strength  in  the  Mohr-Coulomb  model  in  Run  G1  is  the  same  as 
the  constant  shear  strength  in  the  von  Mises  model  in  Run  G2. 
The  Mohr-Coulomb  model  produces  a larger  cavity  (Fig.  3.14), 
more  rapid  attenuation  of  the  peak  stress  for  R > 60  m 
(Figs.  3.16,  3.18,  3.20)  and  greater  radial  displacement 
(Figs.  3,22,  3.24).  For  the  von  Mises  model,  tensile  hoop 
stresses  and  tensile  radial  stresses  occur  (Fig.  3.20). 

Runs  G1  and  G4  both  employ  the  Mohr-Coulomb  model,  but 
the  asymptotic  shear  strength  in  Run  G4  is  only  half  that  in 
Run  Gl.  The  smaller  shear  strength  produces  a larger  cavity 
(Fig.  3.15),  less  rapid  attenuation  for  R > 50  m (Figs.  3.17, 
3.19,  3.21)  and  greater  radial  displacement  (Figs.  3.23,  3.25). 

The  effect  of  introducing  a material  fracture  by  setting 
tensile  stresses  to  zero  that  occur  in  the  von  Mises  model 
is  illustrated  by  comparing  Runs  G2  and  G3.  This  failure 
criterion  causes  a larger  cavity  (Fig.  3.14),  more  rapid  attenua 
tion  of  the  peak  stress  once  the  unloading  wave  from  the  frac- 
tured region  at  the  tail  of  the  pulse  catches  up  with  the 
wave  front  (Figs.  3.16,  3.18,  3.20),  and,  app  arently,  a 
larger  final  radial  displacement  (Figs.  3.22,  3.24).  The 
displacement  histories  also  have  different  shapes  since  the 
imposition  of  the  failure  criterion  inhibits  rebound  (Fig. 

3.24). 
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TABLE  3.5 

GRANITE  PARAMETER  STUDY  RUNS 
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Fig.  3.16--Peak  radial  stress  vs  distance  from  source  for 
Run  G1  (Mohr- Coulomb , high  strength),  Run  G2 
(von  Mises,  without  failure)  and  Run  G3  (von  Mises, 
with  failure) . 
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Fig.  3.17--Peak  radial  stress  vs  distance  from  source  for 
Run  G1  (Mohr- Coulomb , high  strength),  Run  G4 
(Mohr-Coulomb , low  strength),  Run  G5  (capped 
surface,  low  strength),  Run  G6  (Mohr-Coulomb, 
large  cavity)  and  Run  G7  (Mohr-Coulomb,  high 
bulk  modulus) . 
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Mises , with  failure) 


IliL'  tlltii  ol  putting  .-1  cap  on  the  Mohr-Coulomb 
elliptical  yie!,!  surface  ami  introducing  an  associated 
flow  rule  i«  demonstrated  [,»•  comparing  Runs  04  with  Run  (15 

Cap  ",odel  1,1  ,)'1"  (:s  I'roduccs  a slightly  larger  cavitv  ' 

=.140  (ltg.^.15),  much  faster  wave  attenuation  for  R > 30  m 
l-lig-s-  j.l  , a. HI,  3.31),  and  apparently,  smaller  radial 
displacements  (figs.  3.33,  3.35).  The  great  difference  in 
attenuation  ..rise,  from  the  treatment  of  void  collapse  inherent 
n the  cap  plasticity  model,  hut  not  in  the  "oh r- Coulomb 
model  as  treated  in  Run  C4.  The  porosity  effect  could  be 
introduced  m the  p-V-l  equation  ol  state  within  the  'lohr- 
toulomb  model  (or  the  von  Rises  model)  and  a result  closer 
to  that  obtained  by  the  cap  model  would  be  obtained. 

Kun  (it,  uses  the  identical  'lohr-Coulomb  model  used  by 
Ran  (14,  but  the  cavity  in  which  the  S-kT  energy  is  initially 
eposited  ,s  increased  by  an  order  of  magnitude  so  that  the 
initial  pressure  in  the  gas  is  reduced  from  1.54  mbar  to 
hhars.  All  ground  motion  quantities  arc  drastically 
reduced,  at  least  for  1!  < 100  m.  The  gas  pressure  is  too 
near  the  value  of  the  shear  strength  of  the  granite  for  the 
results  to  bo  insensitive  to  the  choice  of  the  cavity  radius. 

Run  (1.  uses  the  same  Mohr-Coulomb  model  as  Run  G4 
except  the  bulk  modulus  is  kept  constant  at  the  value 
derived  from  shock  wave  data  rather  than  permitted  to  vary 
contl nous Jy  with  pressure  according  to  bq . (3.24).  The 
effectively  higher  bulk  modulus  produces  almost  an  identical 

11  lg'  ■5,15)’  bl,t  causes  decreased  wave  attenuation 
and  higher  wave  speed  at  large  distances  (figs.  3.7,  3.19,  3.21). 

'.  ' •"  n>  lt  P roduees  a smaller  radial  displacement 

C-ig.  3.23),  hut  apparently  a larger  radial  d i snlacement  at 
distance  R = 80.54  m from  the  source. 
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3.9  DISCUSSION'  OF  RESULTS 

In  formulating  the  constitutive  equations,  particular 
attention  was  given  to  t he  definition  of  strain  and  strain 
rate.  The  current  definitions  lead  to  constitutive  equations 
consistent  with  those  used  by  Wilkins  and  in  common  use  in  a 
variety  of  spherical,  Lagrangian  codes.  This  is  accomplished 
by  defining  the  principal  strains  as  the  logarithms  of  the 
pi  inciple  stretches,  and  equating  the  strain  rate  tensor 
uith  the  symmetric  part  of  the  velocity  gradient.  Other 
approaches  considered  lead  to  mathematical  and  physical 
difficulties. 

In  particular,  Morland^51^  expresses  the  constitutive 
equation  through  an  equation  of  the  form 

Ur  = al  + bn  + co2  (3.105; 

where  Up  is  the  plastic  part  of  the  right  stretch  tensor. 

If  the  principal  strain  is  defined  as  the  logarithm  of  the 
principal  stretch,  a constitutive  equation  of  this  form  is 
awkward  since  it  contains  the  stretch  explicit'y.  Further- 
more, the  condition  of  incompressibility 

A A A 

1 x 2 . 3 a 

A-  A T = 0 (3.106) 

1 2 3 

becomes  difficult  to  enforce  in  a calculational  scheme.  By 
comparison,  a constitutive  equation  of  the  form 

Dp  = a I + bo 

which  involves  the  stretching,  D,  where 

P = j R (u  U'!  + U-1  fj  | rT 


(3.107) 


(3.108) 


122 


is  very  convenient,  since  the  separation  of  strain  rates 

D = Dp  + De 

follows  trivially  from  the  condition  for  separation  of 
strains,  e . . Furthermore,  the  condition  of 

incompressibility  of  plastic  strains 


aP  aP 

1 A — + — = o 


^ * 
1 2 


aP 

3_ 

aP 

3 


can  be  enforced  by  putting 


(3.110) 


a = -b  tr(c)  , 


(3.111) 


which  makes  the  plastic  strain  rate  proportional  to  the  de- 
viator 4 c stress. 

In  the  approach  taken  by  Clifton  the  principal  elastic 
strains  are  defined  by 


ee  = Ae 
l Ai 


- 1 


(3.112) 


and  the  plastic  strain  bv 


•P  = iP  . 


A i * 1 


(3.113) 


and  the  plastic  and  elastic  stretches  are  related  by 


G p 

Ai  Ai  = > not  summed  on  i . 


(3.114) 


As  in  norland's  approach,  it  is  not  a simple  matter 
to  enforce  plastic  incompressibility  on  the  equations  formu- 
lated in  this  manner,  and  the  constitutive  equations  that 
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would  result  are  quite  complicated  in  their  dependence  on 
stretch  and  strain. 

To  recapitulate,  it  was  found  that  a considerable 
simplification  to  the  constitutive  equations  for  finite 
deformation  plasticity  could  be  obtained  by  defining  the 
principal  strains  as  the  logarithms  of  the  principle 
stretches.  The  strain  rate  tensor  becomes  identical  with 
the  velocity  gradient,  a result  not  always  obtained  with 
previously  suggested  formulations. 

The  constitutive  equation  is  obtained  by  setting  the 
strain  rate  tensor  proportional  to  the  gradient  of  the  plas- 
tic potential.  The  right  stretch  tensor  is  written  as  the 
product  of  an  elastic  and  a plastic  part. 

U = Ue  UP  (3.115) 

and  no  restriction  is  placed  upon  the  plastic  volume.  Conse- 
quently, there  is  a plastic  change  in  volume,  and  it  is  iden- 
tified with  tiie  change  in  pore  volume.  The  change  ir.  volume 
can  be  either  an  increase  or  a decrease,  depending  on  the 
constitutive  equation  and  the  history  of  the  deformation. 

Thus,  the  kinematics  allows  for  either  oulking  or  compaction. 

The  computer  program  and  the  analysis  were  formulated 
in  such  a way  as  to  . llow  a choice  of  either  no  hardening, 
isotropic  hardening  as  expansion  of  the  yield  surface, 
hardening  by  displacement  of  an  elliptical  cap,  or  kinematic 
hardening.  Recent  studies  by  Allen,  et  al.^^  at  S3  and  by 
Maxwell,  et  al_.  J at  Physics  International  have  indicated 
that  the  overshoot  using  the  cap  model  is  less  than  observed 
and  this  conclusion  is  consistent  witli  the  current  results. 

A study  by  McKay  and  Godfrey indicates  that  in  the  "tomb- 
stone" studies  the  attenuation  of  the  calculated  pulse  is  more 
rapid  than  observed.  Both  of  these  deficiencies  in  the 
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theoretical  results  can  be  explained  by  the  absence  of  a Baus- 
chinger  effect  in  the  isotropic  hardening  models.  By  account- 
ing for  the  Bauschinger  effect,  the  velocity  of  rarefaction 
waves  is  reduced  and  the  unloading  takes  place  more  slowly. 

The  velocity  pulse  at  a fixed  station  is  also  stretched  out 
over  a longer  time.  This  is  generally  what  is  required 
to  bring  the  reported  calculations  into  better  agreement 
with  measurements. 


Although  a parameter  study  using  various  plasticity 
models  was  conducted.  Section  3.8,  our  major  emphasis  has 
been  on  trying  to  establish  whether  the  currently  use  1 theories 
are  consistent  with  measurements.  In  view  of  the  disappointing 
correlation  reported  in  the  references  of  the  preceding  para- 
graphs, it  seemed  worthwhile  to  make  an  overall  examination 
of  the  existing  data  and  the  material  model.  It  was  con- 
cluded that  permanent  displacement  seems  to  be  the  most 
reliable  measurement,  and  good  correlations  were  generally 
found  between  measurements  and  the  simple  formula,  6 = R3/3R2. 
One  Piledriver  displacement  is  a notable  exception.  This 
displacement  lies  a factor  of  three  above  the  "best  fit"  by 
Borg.  Its  reliability  is,  however,  questionable  simply  in 
view  of  the  substantial  difference  between  the  velocity  and 
integrated  acceleration  measurements  at  the  1543-ft  station, 
reported  by  Ferret.  (Though  the  two  measurements  agree  quite 
well  at  early  times,  they  diverge  after  0.3  sec  real  time, 
or  76  msec  in  the  scaled  plot  of  Fig.  3.10.  The  accelerometer 
trace  would  lead  to  final  values  of  displacement  well  below 
the  best  estimates  given  by  Borg,  and,  in  fact,  would  ulti- 
mately have  the  wrong  sign  since  the  duration  of  the  negative 
velocity  phase  is  indefinitely  long.) 

The  peak  velocity  and  final  displacement  predicted  by 
the  "cap"  model  are  too  low.  This  can  be  interpreted  as  a 
consequence  of  the  excessively  high  rate  of  unloading  when 
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the  Bauschinger  effect  is  not  accounted  for.  A more  realis- 
tic attenuation  is  predicted  with  the  kinem.tic  hardening 
mode  1 . 

There  remain  a number  of  specific  investigations  which 
should  be  carried  through  to  explain  the  existing  discrepancies 
between  calculation  and  measurement.  Of  particular  concern 
is  the  reduction  in  flow  stress  by  a factor  of  six  that  is 
lequired  to  bring  calculations  into  agreement  with  shot  data. 
The  existence  of  large  cracks  and  pore  water  are  probably 
responsible  for  this  discrepancy.  It  would  be  desirable  to 
get  quantitative  information  on  this  question,  and  on  the 
relative  importance  of  pore  water  and  size  effect  on  strength. 

The  effect  of  material  dilatancy  has  not  been  examined 
in  enough  detail.  Dilatancy  was  considered  at  one  time  to  be 
one  of  the  main  effects  causing  discrepancy  between  calculation 
and  shot  data.  The  effect  may  be  important,  but  as  a result 
of  the  current  study,  it  is  believed  that  dilatancy  should 
be  considered  in  conjunction  with  a model  that  accounts 
for  the  Bauschinger  effect. 

The  lithostatic  pressure  in  the  far-field  causes  the 
cavity  to  rebound.  In  addition,  it  significantly  influences 
the  strength.  We  have  not  accounted  for  either  of  these 
depth  effects  in  the  current  calculation. 

Finally,  we  note  that  a reasonable  model  for  gross 
fracture  due  to  excessive  strain  or  tensile' stress  should  be 
added  to  the  simulation.  The  treatment  of  fracture  should 
reflect  the  finite  time  required  for  cracks  to  propagate,  and 
for  the  material  to  lose  its  competence.  Such  a model  has 
not  yet  been  incorporated  into  the  computer  programs  commonly 
used  to  calculate  underground  nuclear  explosions,  but  it  is 
felt  that  the  delay  in  fracturing  and  the  tesidual  strength 
would  significantly  influence  the  late  behavior.  Without  a 
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gradual  failure  model,  the  material  would  tend  to  form  iso- 
lated shells  in  the  calculation,  whereas  in  practice  some 
competence  is  retained  almost  everywhere. 


IV.  TllliORY  OF  1 NTHRACTING  CONTINUA 


4.1  I NT kODUCT I ON 

Shock  wave  propagation  in  composite  materials  (e.g., 
fiber  reinforced  composites,  water  saturated  rocks)  is  of 
interest  to  many  areas  of  science  and  engineering.  Composite 
materials  display  certain  effects  (e.g.,  geometric  wave  dis- 
persion, internal  dissipation)  which  cannot  be  adequately 
modeled  within  the  usual  restrictions  of  a homogeneous,  iso- 
tropic media.  The  material  properties  of  the  various  consti- 
tuents, geometrical  arrangement  and  the  porosity  (or  void 
content)  all  affect  the  thermo -mechanical  response  of  the 
composite.  Thus,  it  is  necessary  to  develop  an  analytical 
model,  based  on  the  microstructure  of  the  composite  which 
will  enable  one  to  evaluate  the  thermo -mechanical  behavior  of 
the  material  under  various  static  and  dynamic  loading  condi- 
tions. However,  the  variability  in  material  properties  and 
the  difficulty  of  characterizing  interfaces  make  a purely 
microscopic  approach  somewhat  unattractive.  A practical 
ana/  ^ical  model  should  provide  an  average  description  of 
the  constituents  rather  than  a detailed  thermo -mechanical 
description  at  each  material  interface  at  each  instant  of 
time.  The  theory  of  interacting  continua  (TINC)  provides 
a means  for  proceeding  directly  to  the  desired  macroscopic 
level . 

In  3SR-267^^  and  3SR-648,^^,  this  theory  was 
introduced  to  provide  a framework  for  describing  the  behavior 
of  a geologic  (dry  and  partially  saturated  tuff)  mateiial  in 
terms  of  the  isolated  behavior  of  the  constituents.  The  major 
emphasis  therein  was  placed  on  developing  a mechanical  theory 
appropriate  for  planar  stress  wave  propagation.  This  model 
was  incorporated  into  the  planar  POROUS  code.  Additionally, 
Hugoniot  relations  were  derived  for  a binary  mixture  in 

Preceding  page  blank 
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3SR-648.  During  the  past  year,  the  prime  objective  has 
been  to  develop  a thermodynamical  theory  suitable  for  wave 
propagation  studies  in  planar  or  spherical  geometry.  After 
outlining  the  conservation  relations  in  Section  4.2,  we  dis- 
cuss the  various  interaction  terms  in  Section  4.3.  The  con- 
stitutive relations  are  given  in  Section  4.4.  The  incorporation 
of  the  thermodynamical  theory  into  a completely  new  version  of 
POROUS  treating  both  1-D  planar  and  spherical  configurations 
is  discussed  in  Section  4.5.  Finally,  in  Section  4.6,  we 
discuss  some  material  parameter  calculations  using  this  code. 
Only  limited  calculations  have  been  made  since  the  new  POROUS 
code  is  still  in  its  final  development  stage  at  the  time  this 
report  is  being  written. 
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4,2  fflNHRAL  CONSI-RVATlnw  uwc 

composite  body ’is  occupieTby  particle"*  f"311  VOl™e  °f  'he 
A (c  ■ 1,  2,  ...  r,  ,,  7 Particles  of  each  constituent 

has  a velocity  field  ’ ;x  "T0  ’ constituent  (?J 

i denotes  the  positio„“Vecior  i„ 'T  ^ C°mp°Site  Khe™ 
fixed  Newtonian  frame)  and  t i .(Wlth  reference  to  a 

constituent  per  . enotes  time.  The  mass  of 

partial  density  (p3  rx  V°  Un*  °f  co^POsite  is  called  its 

volume  of  composit  ™ ^ Per  unit 

PU,tj  is  given  by: 


imilarl) , the  total  stress  tensor  a 
area  of  the  composite  can  be  decomposed  S"°Clated . Wi th  a unit 
e associated  with  each  component  ^ P*rtl91  StrCSSes 


(For  a physical  explanation  of  n,„ 

herein,  reference  is  made  to  a paper^^T 

momentum^nd^nergy 'ma^be^ritte^for11^6'!!31^^  °f  ’ 

and  3SR-648J  . „ will  „ow  assume  that ^ ! f3SR-267 

forces  are  absent,  (3)  there  is  no  heat  lf  ^ 

external  world  m nn  * transfer  to  the 

constituents  due  to  chelical'j'nt^"  b°tKeen  thc' 

stress  tensors  are  symmetrT^Vith  T ^ 

the  conservation  relations  become:  S aSsumPtions. 
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Cont inuity : 


a 


(nH“)  M (a) 

- ut'p— ■ + P div  v 


dV 


/ 


,(“J  r MM 

-sf-  dv  * / 0 vj  nj  ds 

v s 


= 0 


Momentum : 


/ 


(a)  (a) 

(a)  D v- 
p — Bt-1  dv 


r , / (a)  (a)  v /•  (a)  (.aj  (.aj 

' / ! T ( " vi)dV  " / d vi  vj  "j 

*/V  ‘'s 


(a)  (a)  (a) 


/(a)  f (a) 

a n^  dS  + I p dV 


Energy : 

/ 


w <;>  ' 

p T5T 


(a)  i (a)  (a) 

E + 7 vj  Vj 
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'(a). 

((a) 

i (a)  (a)  \ "] 
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Jv 

P 1 

i E 

VJ 

dV 

r (a) 

f(ot) 

1 

(a)  (a)  “I 

(a) 

* f 0 

LE 

+ 7 

vj  vj  J 

Vi 

n. 

i 

dS 


(4.3) 


dS 


(4.4) 
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r f 


/(“)  O)  /•  (a) 

°ij  vi  nj  ds  - J q-i  nj 


dS 


f ( (a)  (a)  (a) ) 

I )0  0 • V + p ip  dV 


(4.5) 


where 


(a) 

v 

(a) 

P 

(a) 

aij 

(a) 

P 6 


(a)  (a) 

D / Dt  = 9/3t  + v • grad 

velocity  of  material  ^ 
partial  density  of  ^ 


(a) 


(a) 

E 

(a) 
P p 


n 


(a) 

3. 


= partial  stress  tensor  for  a 

= momentum  supply  to  per  unit  mass 

of  composite  due  to  interaction  forces 

• ■ . fa) 

specific  internal  energy  of  4 

(a) 

= energy  supply  to  4 per  unit  mass  of 

composite  due  to  interaction  with  the  remaining 
a-1  constituents 

= normal  to  the  surface  S 

heat  flux  vector  into  ^ from  the  remaining 
a-1  constituents 
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The  requirements  that  the  total  momentum  and  energy 
contributions  of  the  internal  material  interaction  forces  be 
zero  may  be  written  as 


a=  1 


(4.6) 


'/(a)  (a) 

(a)\ 

(«)' 

\ • v + 

* ) 

- div 

= 0 

(4.7) 

a= 


The  balance  equation  for  internal  energy, 
tained  by  combining  Eqs . (4.3),  (4.4),  and  (4.5). 


(a) 

E , may  be  ob- 


3 

3t 


(a)  (a)  (a) 
P E vi 


ni 


dS 


dS 


+ 


(a) 

p ^ dV  . 


(4.8) 


In  writing  down  the  above  conservation  relations  in  the 

TINC  framework,  no  reference  is  made  to  the  actual  mean  area, 

(a)  (a)  (a) 

m , and  mean  volume,  n , occupied  by  <5  per  unit  cross- 
sectional  area  and  volume  of  the  composite.  If  the  inter- 
action terms  and  constitutive  relations  for  the  composite  are 
to  be  expressed  in  terms  of  the  behavior  of  the  isolated  consti- 
tuents, reference  must  be  made  to  the  actual  constituents.  In 

(a)  e 

3SR-267,  effective  densities,  p , and  effective  stress 
(a)  . 

tensors,  a e,  were  defined  in  terms  of  partial  densities 
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and  partial  stress  tensors  by  the  scaled  relations: 


(a)  (a)  (a) 

p = n p 

(a)  (a)  Cot) e 

o = m a 

whe  re 


a= 1 a=l 


(4.9a) 

(4.9b) 

1 (4.9c) 


If  the  composite  is  isotropic  so  that  each  plane  through  the 

(a) 

medium  intersects  the  same  area  fraction  of  6 . the  area 

(a) 

and  the  volume  fractions  are  the  same  for  each  & . 


0)  (a) 

m = n (4.10) 


In  this  case  a single  scaling  function  occurs  in  the  relations 
(4.9). 

In  the  following  discussion  a = 1,  2,  3 will  be  used 
to  designate  rock  (poreless),  water  and  voids,  respectively. 
The  interaction  *:erms  in  the  next  section  are  derived  for  the 
case  when  n = 0,  i.e.,  a fully  saturated  porous  solid. 

The  applicability  of  these  terms  to  the  unsaturated  case  is 
discussed  in  Section  4.4. 
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4.3  INTERACTION  TERMS 


Tor  sake  of  convenience,  we  will  restrict  the  discussion 
in  this  section  to  plane  wave  propagation  in  the  x-direction. 
The  interaction  terms  derived  herein  are,  however,  also  appli- 
cable to  spherical  wave  propagation.  For  planar  wave  propa- 
gation, the  conservation  relations  (4.3),  (4.4)  and  (4.8)  may 
be  written  in  the  differential  form  as: 


(a) 


(°)  (a)  (c° 


P 1 


d P ♦ V 

3 1 

d p 

3x 

«J  v _ n 

■ * 0 3x  - ° 

(4.11) 

(a)  (a) 

v ♦ v 
3 1 

(a)  v 

3 V 
3x  i 

(a) 

1 (°0  9 
I-  <=  6 * -53T 

(4.12) 

^ («) 
TT*  v 

(a), 
3 E ) 
3x  j 

(a)  («) 

1 * °x  -75T  * p * 

(4.33) 

In  writing  Eq . (4.13),  we  assumed  no  heat  transfer  between  con- 
stituents , 

(a) 

q = 0 (4.14) 


Our  primary  interest  is  in  mixtures  in  which  the  second 
component  is  a fluid. 


(2)  (2) 

a..  = -p  6 


(2)  (2) (2) 

P = n p 

Rewriting  Eq.  (4.6),  we  obtain 

(1)  (2) 

p 3 = -p  3 = p3 


(4.15a) 

(4.15b) 


(4.16) 
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The  interaction  term,  p$,  contains  both  dilatational  and 
diffusive  parts 


pe  - P0d  + pQn  (4.17) 

To  derive  an  expression  for  p£3,  we  consider  the  ordinary 
momentum  equation  for  fluid  flow  with  friction  through  a variable 
area  tube , 


(2) (2) 


m 


/ (2) 
a v 
\~TT 


(2) 

v 


(2) 
a v 
"3x” 


pon 


Cot) 

m 


(2) 


Tx 


(4.18) 


where  -p^n  denotes  the  drag  experienced  by  the  fluid  as  it 
flows  through  the  rock.  Combining  Eqs . (4.9),  (4.10),  (4.12), 
(4.10),  (4.17)  and  (4.18),  we  obtain: 


pe 


(2) 
a n 
3x 


♦ P n 

0 


(4.19) 


Note  that  Eq . (4.19)  holds  only  as  long  as  (4.10)  is  true. 

We  will  now  assume  that  the  diffusive  force  p n depends 
on  the  two  velocity  fields  in  the  following  manner: 

I (2)  (1)\ 

Pon  = Pod  \ v ' V / C4>2°) 


where  d has  the  dimension  of  the  reciprocal  of  time.  Note 
that  d is  rela.ed  to  Darcy's  law  (3SR-648) 


P d 

o 


(4.21) 


where  y and  k denote  the  kinematic  viscosity  of  the  fluid 
and  the  permeability  of  the  rock,  respectively.  For  water,  y 
is  equal  to  0.01002  g/sec-cm.  For  tuff,  k is  in  the  range 
50-10,000  ydarcies  (1  darcy  = 0.987  10‘8  cm2). 
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It  appears  worthwhile  to  note  here  that  (1/d)  is  a 
measure  of  the  momentum  relaxation  time,  i.e.,  the  time  re- 
quired to  exchange  momentum  between  the  two  materials.  In 
numerical  computations , this  fact  can  be  very  important  in 
deciding  the  allowable  time  steps.  Ideally,  the  time  step 
should  be  an  order  of  magnitude  less  than  the  momentum  lOlaxa- 
tion  time,  i/d. 

(a) 

4*3.1  Interaction  f.nergy  Term:  p \p 

I Provided  there  is  no  heat  flow  between  the  constituents 

1 q = 0 ),  the  requirement  that  the  energy  contributions  of 
the  interaction  forces  be  zero  may  be  written  as 


/Cl)  (2)  \ /Cl)  C2) v 

P3  \ v - v p ( \p  + ip  j = o (4.22) 

(2) 

The  internal  energy  balance  relation  for  6 can  be 
rewritten  as: 


p Dt~ 


(2) 

E 


(2) 

-P 


(2) 
3 v 
9x 


(2) 

+ p ip 


(4.23) 


noninteraction  energy  ^rm,  p ip  , contains  both  dilatational , 
o j , and  diffusive,  p \pg  , contributions 


(2)  (2)  (2) 

p ip  = p + P iPs  (4.24) 


Also,  we  have 


(2).  (2) 

E e s E 


To  derive  an  expression  for 


(4.25) 

we  consider  the  usual  internal 


k. 
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balance  equation  for  fluid  flow  with  friction 

(25 (2)  C2J  (2)  (2) (2)  (2),  (2) 

n p E = -n  p div  v ^ + p ips  (4.26) 

(2) 

where  p ip  denotes  the  diffusive  energy  contribution  to  the 
fluid  as  it  passes  through  the  rock. 

Now,  we  have  from  mass  conservation: 
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div  v e 


C2)(2) 

-Dp' 

ITT- 


/u, 

/ P 


(2) 


Dt 


/ P 


+ m 


n 


(2) (2) 
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T77 
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C2)  (2) (2) 

9 v 1 D n 

9x  (2) 

n 


Tt 


(1.27) 


Substituting  from  Eq . (4.27)  into  Eq.  (4.26)  and  utilizing 
Eqs.  (4.22)  and  (4.24),  we  obtain 


(2) 
P *d 


(2) 
9 n 
9 1 


(2) 

v 


(2) 

9 n 

9x 


(4.28) 


The  requirement  that  the  energy  contribution  of  the 
internal  interaction  forces  be  zero  yields 


(2>, 

- P 


(2) 

9 n 
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(2)  Q-) 

9 n UJ 

'Tjr  v 


(i) 

p 'P 


Pon 


/CD  (2)v 

( v - v ) 


(2) 

P 


= 0 (4.29) 
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can  be  split  up  into 


The  internal  energy 

(] ) 

two  parts,  p \p. 


contribution 
. (1) 
and  p \ps  . 


(1) 

P ij> 


(i)  (1)  C1) 

p -tj  = p -^d  + p \ps  (4.30) 


We  have  thus  twcjo  ^qua^iyns , £4 . 29)  and  (4.30),  with  four  un- 
knowns, p ^ , p ip^,  P p i's.  It  will  now  be  assumed  that 

the  dilatation  and  diffusive  contributions  must  be  separately 
equal  to  zero , i . e . , 


CD 
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CD, 


C2) 
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(4.31) 


CD  C2)  / (2)  (1)  v 

P + p = pQn\  v ' v ) (4.32) 


and  incompressible  materials  ( n = constant).  However,  in  the 
general  case,  it  has  to  regarded  as  a constitutive  assumption. 

We  still  requii^  another  constitutive  assumption  to 
separately  evaluate  \ps  and  $ . Since  diffusion  is  a 

dissipative  process,  the  following  inequalities  follow: 


CD  (2) 

P i|»s  > 0 P > o 


We  now  assume  that 


is  zero. 


Then 


(2) 
P *s 


p n 

0 


(2)  (ID 

v • v > 


(4.33) 


(4.34) 


The  last  assumption  states  that  the  fluid  receives  all  the 
diffusive  energy  contribution,  and  is  justified  by  the  fact 
that  the  thermal  conductivity  of  the  rock  is  usually  much  less 
than  that  of  the  fluid  [see  also  Ref.  54]. 


A more  general  procedure  for  evaluating  the  various 
interaction  terms  is  outlined  in  the  Appendix. 

4.3.2  Relation  to  Earl * . Vork 

In  previous  work  (3SR-648)  , it  was  assumed  that 


Cot) 
p Ip 


Cot) 
■P0  v 


(4.35) 


If  one  now  assumes  that  the  present  splitting  of  p into 

dilatation  and  diffusive  components  is  valid,  then  it  follows 
from  Eq.  (4.35) 
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(4.36) 


(4.37) 


It  is  easily  seem  that  Eqs . (4.36)  and  (4.37)  cannot  possibly 
satisfy  inequalities  (4.33).  Hence,  it  is  concluded  that 
Eqs.  (4.36)  and  (4.37)  are  invalid.  It  may  also  be  verified 
that  for  the  steady  case,  present  expressions  for  p CiR  and 
P Reduce  to  the  Eq.  (4.35).  The  present  expressions  for  p 
fif  _p(fj  are  identical  with  Eqs.  (4.36)  and  (4.37)  for  S 

v » p0n  - 0 and  v =0.  Thus,  we  observe  that  the 
Hugoni^anajjjis  of  3SR-648  remains  valid  except  in  the  case 
when  v f v and  P(jn  ^ 0.  The  results  of  3SR-648  may  be 
modified  in  a straightforward  manner  to  include  this  case  as 

well.  We  will  not,  however,  pursue  this  question  any  further 
here . 

The  derivation  of  interaction  terms  for  spherical  flew 
is  similar  to  the  planar  case.  In  fact,  Eqs.  (4.19),  (4.28) 
and  (4.34)  are  directly  applicable  to  the  spherical  case  with 
x replaced  by  the  spherical  space  vaiiable  r. 
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4.4  CONSTITUTIVE  RELATIONS 


In  3SR-648  mechanical  constitutive  relations  were 
presented  for  a partially  saturated  rock  under  the  assumption 
of  disconnected  voids  (voids  contained  in  rock  matrix  and 
separated  from  saturated  pores).  The  major  effort  this  last 
year  has  been  directed  towards  (1)  improving  the  deviatoric 
stress -strain  relations  for  the  rock,  (2)  including  thermal 
effects  into  constitutive  laws,  and  (3)  developing  a suitable 
thermal  crushup  model.  This  work  is  discussed  more  fully  in 
the  following  subsections. 


4.4.1  Improved  Formulation  for  Deviatoric  Stresses 


In  the  previous  TINC  work,  the  following 
relations  were  employed  for  the  rock  component 
elastic  regime: 


constitutive 

(D  • 

4 in  the 


a)  a)  a) 

°i  ■ • p * si 


(4.38a) 


(4.38b) 


(1) 
= 2 n 


■ft 

' o 


(4.38c) 


(1) 


Here  'o^,  ^ and  '‘S^  denote,  respectively,  the  partial 
stress,  the  partial  pressure  and  the  partial  deviatoric 
stresses . 


% 


is 


A^'s  are  the  partial  principal  stretches;  ) 

the  partial  Jacobian  of  deformation;  n is  the  volume  frac- 
(1) 

tion  of  4 . Functional  P determines  the  pressure  response 


of  the  pore  less  rock;  p 
less  rock. 


is  the  shear  modulus  of  the  pore- 


For  a dry  rock,  partial  quantities  correspond  to  the 
quantities  actuallj^measured  in  a laboratory.  Thus,  for 
example,  X and  X^  are  related  to  the  bulk  volumetric 
strain,  G,  and  the  bulk  principal  strains  c^'s  through 
the  relations  : 

x = i + e 

X.  = 1 + e.  (4.39) 

G = c . + c + c 

1 2 3 

We  now  assume  th^t  for  the  dry  rock  in  the  regime  of  small 
deformations,  n"  depends  only  upon  the  volumetric  strain, 

(D  fl)  / , v 

n = n^  |l  + a^G  + a Gz  + ...  J (4.40) 

where  a , a , ...  are  constants.  Substituting  from  Eqs . (4.39) 
1 2 

and  (4.40)  into  Eq  . (4.38)  and  neglecting  terms  of  0(61 2),  there 
follows : 

(1)  (1) 

p = - n K (1  + a ) * - KG  (4.41) 

o s 1 

(1)  (1)  (1) 

S.  - S,  = 2 n p (e.-e,)  = 2p(e. - e - ) (4.42) 

1 J o i J-  J pij 

where 

K = bulk  modulus  of  poreless  rock 
K = bulk  modulus  of  porous  rock 
Pp  = shear  modulus  of  porous  rock 
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Equations  (4.41)  and  (4.42)  imply  that: 


(1) 

K - n K (1  + a ) 
o s 1 

(1) 

„ = n u. 


(4.43) 

(4.44) 


It  is  shown  elsewhere  ^55  ^ that  relationship  (4.43)  is 
consistent  with  the  linear  theory  of  elasticity.  However,  we 
note  here  that  Eq . (4.44)  represents  a physically  unrealistic 
result.  Most  porous  rocks  exhibit  a very  small  shear  modulus 
at  low  confining  pressures  . Furthermore,  the  measured  values 

of  the  shear  modulus  have  been  found  to  be  path-dependent . ^57  »58  ^ 

The  relationship  (4.44)  is  a consequence  of  the  assumed 
relationship  between  the  effective  deformation  gradient  p e 

• \ CX  J ~ 

and  partial  deformation  gradient  Y-  . 


and 


(4.46) 


Equation  (4.46)  is  identically  true.  The  justification  for  Eq . 
(4.45)  is  not  so  clear.  There  is  little  reason  to  expect 
that  A^'s  should  be  related  to  A^'s  through  an  isotropic 
relation  of  the  form  of  Eq.  (4.45).  As  an  example,  let  us 
consider  a bar  with  a spherical  cavity  subjected  to  uniaxial 
loading.  The  spherical  cavity  would  deform  into  an  ellip- 
soidal cavity.  However,  for  Eq.  (4.45)  to  apply  the  spherical 
cavity  must  retain  its  shape.  We  remark  here  that  Eq.  (4.45) 
is  approximately  correct  for  water  saturated  low  strength 
rocks  like  tuff.  In  low  strength,  high  porosity  rocks  the 
pore  pressure  is  nearly  equal  fin  confined  tests)  to  the 
confining  pressure  and  helps  to  retain  the  original  pore 
shape.  However,  for  dry  rocks  and  high  strength  saturated 
rocks,  Eq.  (4.45)  can  lead  to  gross  error^.^  Instead  we  intro- 
duce the  following  relationship  between  F e and  ^ : 


Ae 


(4.47) 
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where 


Z + Z + Z s 1 

1 2 3 


(4.48) 


The  restriction,  Eq.  (4.48),  on  the  exponents  Z^  is  required 
to  satisfy  Eq . (4.46).  With  the  decomposition  relation  in 
Eq.  (4.47),  the  deviatoric  stress  strain  relationship,  replacing 
Eq . (4.38),  becomes  : 


where  the  repeated  subscripts  i and  j do  not  imply 
summation . 


Substituting  from  Eqs.  (4.39)  and  (4.40)  into  Eqs . (4.49), 
we  are  led  to  the  following  relationship  between  y and  y : 


(1) 


n y j 
o i ) 


Ei'Ej 


iEi 


■■i  * V(£rcPj 


(4.50) 


Equation  (4.50)  implies  a path-dependent  relationship  between 
u j and  Up.  Thus,  for  example,  for  uniaxial  strain  and  tri- 
axial  tests  we  have: 


Uniaxial  Strain: 

e i 0 , c = c = 0,  0 = c , 

1 2 3 1 

' £,  ' C1  - V/2 

(l)  ( | 

up  = \ M1  * a,(-1/2  * 3V2)| 
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Triaxial  Loading: 


e^O,  £=£^0,  6 = e + 2e 

12  3 12 

£,)/2 

(D  I E ♦ 2e  I 

"p  * \ ",  lta,  e - c 2 t'^2  * 3V2> 

The  factors  will  in  general  depend  upon  the  deformation 

path  and  may  be  evaluated  if  appropriate  data  are  available. 
This  could  be  quite  a cumbersome  process. 

Fortunately,  for  geologic  composites,  only  the  rock 
matrix  can  sustain  shear  stresses,  and  the  partial  shear 
stresses  for  the  rock  equal  the  composite  shear  stresses. 

It  is,  therefore,  unnecessary  to  calculate  effective  shear 
stresses.  Instead,  one  can  directly  postulate  a shear  law 
of  the  form: 

CD  CD 

s.  - s.  * 2U  (X.  - Xj)  (4.51) 

As  was  remarked  earlier,  up  will  in  general  depend  upon  the 
deformation  path  and  pore  pressure.  We  shall  presume  here 
that  such  a functional  dependence  is  known  from  the  experi- 
mental data.  Thus,  the  mechanical  TINC  constitutive  relations 
become  Eqs . (4.38a),  (4.38b),  and  (4.51)-  The  new  formulation 
retains  the  major  advantage  of  the  previous  work,  viz.  the 
treatment  of  porosity  as  an  independent  variable.  This  is 
quite  important  in  studying  the  effect  of  pore  pressure. 

The  constitutive  law  for  shear  stresses,  Eq . (4.51),  sidesteps 
the  difficulties  inherent  in  relating  the  components  of  the 
partial  deformation  gradient  to  those  of  the  effective  defor- 
mation gradient.  It  is  also  to  be  noted  that  Eq.  (4.51)  uti- 
lizes the  information  which  can  be  readily  obtained  in  the 
laboratory . 


4.4.2  Thermodynamic  Constitutive  Laws 

For  a thermodynamic  theory,  we  have  to  develop  consti- 
tutive laws  for  rock  partial  stresses,  ^a],  and  fluid  partial 
pressure,  p . As  n^J^d  earlier,  may  Le  decomposed  into 

an  isotropic  part,  - p , and  a deviatoric  stress,  sj>  Eq.  (4.38). 
We  will  assume  that  Eq.  (4. SI)  for  shear  holds  for  the  thermal 
case  as  well.  If  nectary,  the  shear  modulus,  p , can  be 
made  to  depend  upon  E and  ^E^  as  well.  For  plastic  flow, 
we  will  assume  von-Mises'  law 
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(l)2 
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(4.52) 


where  denotes  the  yield  stress  of  the  porous  rock  in 

simple  tension  or  compression.  The  yield  stress,  Y , may 
depend  upon  the  two  pressures,  p and  , and  the  inter- 

nal energies,  E^  and  . We  will  not,  however,  attempt 
here  to  specify  any  particular  functional  relationships  for 

MP  3nd  YP’  but  wil1  assume  these  to  be  known  from  the 
experimental  data. 

( a) 

The  pressure,  pa,  for  4 constituent  in  isolation 
is  a function  of  density  and  internal  energy. 


P0  * VP.  V E> 


(4.53) 


Alternately,  p^  may  be  regarded  as  a function  of  p,  p and 
T (temperature).  . Specific  functional  forms  for  NTS  tuff  and 
water  are  discussed  elsewrhere  in  this  report. 

Within,  the  mixture,  the  isotropic  part  of  the  stress 
(a) 

tensor,  - p , is,  therefore,  given  by 
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This  completes  our  discussion  of  the  thermodynamics  except 
insofar  as  we  still  need  to  outline  a procedure  for  determining 
n . This  is  considered  in  the  next  section. 


4.4.3  Crushup  Model 


In  3SR-648,  mechanical  crushup  relations  were  presented 
for  dry  and  partially  saturated  tuff.  During  the  present 
contract  period,  simpler  crushup  relations  have  been  developed. 
These  are  discussed  in  detail  in  Section  2.2.3  of  this  report. 
However,  for  the  sake  of  completeness,  we  will  outline  the 
model  here  in  TINC  notation. 

(l)/(2)FOr  thf3?ry  pof?ys  rock’  we  have  onlX  one  unknown, 
n \ n - 0 > n = 1-  n j.  It  is  useful  to  define  a new 
variable,  p such  that 


/d)\  1 

p / " TU 


A quadratic  relationship  is  assumed  for  ct(  p ) : 

/d)\  / (1)  \2 

a(  p ) - 1 + ( a a ' 1 ) ^ 1 ' P /Pc  j 


(4.55) 


(1)  3(d5 

for  p < pc  and  2L£_  > 0 

/ d)  \ (1) 

«\  P I = 1 for  P > pc  (4.56) 

On  unloading,  a is  held^constant . Here  denotes  the 

initial  value  of  a|=  1/  n^  j and  p^  is  the  crushup  pressure. 
Thus,  the  only  required  inputs  for  this  model  are  a and  p . 


149 


For  the  par^j.^lly  rocj^  has  two  unknown 

volume  fractions,  n and  n ( n = 1-  n - n ).  In  3SR-648,  two 
postulates  concerning  pore  crushup  were  introduced.  The  pores 
may  be  considered  to  be  connected  such  that  the  water  can  move 
freely  between  the  pores.  For  shock  wave  studied,  however,  the 
more  appropriate  postulate  is  that  of  disconnected  pores  that 
are  either  completely  water  saturated  or  gas  filled  with  no 
water.  In  this  model  the  partially  saturated  rock  may  be  con- 
sidered to  be  a composite  in  which  the  components  are  water  and 
distended  rock. 

The  disconnected-pores  hypothesis  will  be  used  here  to 
develop  the  crushup  relations  (see  also  3SR-648).  We  now  de- 
fine a to  be: 


(1)  (3)  (2) 

_ n + n 1 - n 

' ~cn ut~ 

n n 


(4.57) 


The  crushup  relation  (4.56)  is  modified  as  follows: 

a ~ 1 + (a/1)  (1-P/PC)2  (4.58) 

where 
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CD 
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~nr 
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(4.59) 


We  need  an  additional  relation  to  complete  the  algebraic  loop. 
This  is  obtained  by  equating  the  effective  pressures  in  the 
water  and  the  distended  tuff  components,  i.e., 
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This  completes  the  description  of  the  crushup  model. 


L 
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It  is  now  readily  seen  that  for  the  disconnected 
pores  postulate,  the  interaction  terms  derived  in  Section 
4.3  for  the  completely  saturated  rock  are  directly  appli- 
cable to  the  partially  saturated  case. 
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4.5  THERMODYNAMIC  POROUS  CODE 

A description  of  the  mechanical  POROUS  code  finite 
difference  scheme  has  been  given  in  3SR-648;  this  formu- 
lation did  not  include  thermodynamics,  i.e.,  the  effects  of 
internal  energy  on  the  equations  of  state  of  the  constituents. 
The  finite  difference  scheme  used  in  the  non-thermodynamic 
POROUS  code  was  based  on  the  Lax -Wendrof f technique,  and 
furthermore  was  not  written  in  conservation-law  foim.  The 
thermodynamic  version  of  POROUS  incorporates  the  effect  of 
internal  energy  on  pressure  in  the  two  constituents,  uses  a 
finite  difference  scheme  based  on  the  so-called  "leap-frog" 
technique  and  is  written  in  conservation- law  form,  so  that 
the  masses  of  the  constituents  and  total  energy  are  exactly 
conserved  (except  for  roundoff  error),  and,  in  planar  geometry, 
momentum  is  also  conserved.  It  is  well  known  that  the  Lax- 
IVendrof f scheme  often  results  in  "overshoots"  of  dependent 
variables,  such  as  pressure  and  density,  at  a shock  front;  the 
use  of  the  "leap-frog"  scheme  tends  to  eliminate  these  over- 
shoots, which  can  be  particularly  troublesome  in  problems 
involving  shocks  at  low  stress  levels  (several  kbar  or  below). 


In  the  POROUS  treatment  of  water-tuff  mixtures,  only 
the  tuff  has  material  strength;  hence,  it  is  convenient  to 
employ  the  conservation  laws  for  mass,  momentum  and  energy  in  a 
coordinate  system  w'hich  moves  with  the  tuff,  i.e.,  coordinates 
which  are  Lagrangian  with  respect  to  the  tuff.  In  general, 
there  will  be  transport  of  water  from  cell  to  cell,  but  the 
mass  of  tuff  in  a cell  remains  constant.  In  this  coordinate 
system,  the  conservation  laws  of  Section  4.2  appear  as  follows: 
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4.5  THERMODYNAMIC  POROUS  CODE 

A description  of  the  mechanical  POROUS  code  finite 
difference  scheme  has  been  given  in  3SR-648;  this  formu- 
lation did  not  include  thermodynamics , i.e.,  the  effects  of 
internal  energy  on  the  equations  of  state  of  the  constituents. 
The  finite  difference  scheme  used  in  the  non- thermodynamic 
POROUS  code  was  based  on  the  Lax --Wendrof f technique,  and 
furthermore  was  not  written  in  conservation - law  form.  The 
thermodynamic  version  of  POROUS  incorporates  the  effect  of 
internal  energy  on  pressure  in  the  two  constituents,  uses  a 
finite  difference  scheme  based  on  the  so-called  "leap-frog" 
technique  and  is  written  in  conservation- law  form,  so  that 
the  masses  of  the  constituents  and  total  energy  are  exactly 
conserved  (except  for  roundoff  error),  and,  in  planar  geometry, 
momentum  is  also  conserved.  It  is  well  known  that  the  Lax- 
Wendroff  scheme  often  results  in  "overshoots"  of  dependent 
variables,  such  as  pressure  and  density,  at  a shock  front;  the 
use  of  the  "leap-frog"  scheme  tends  to  eliminate  these  over- 
shoots, which  can  be  particularly  troublesome  in  problems 
involving  shocks  at  low  stress  levels  (several  kLar  or  below) . 

4.5.1  Conservation  Law  Equations  in  Moving  Coordinates 

In  the  POROUS  treatment  of  water-tuff  mixtures,  only 
the  tuff  has  material  strength;  hence,  it  is  convenient  to 
employ  the  conservation  laws  for  mass,  moment'.im  and  energy  in  a 
coordinate  system  which  moves  with  the  tuff,  i.e.  , coordinates 
which  are  Lagrangian  with  respect  to  the  tuff.  In  general, 
there  will  be  transport  of  water  from  cell  to  cell,  but  the 
mass  of  tuff  in  a cell  remains  constant.  In  this  coordinate 
system,  the  conservation  laws  of  Section  4.2  appear  as  follows: 
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Continuity : 
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Energy : 
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The  ordinary  time  derivative,  d?,  is  used  in  the  above  equations 
because  one  is  evaluating  the  time  rate  of  change  of  mass,  momen- 
tum and  energy  in  a cell  of  nonzero  width;  the  hydrodynamic 
derivative,  ^ is  appropriate  to  locally  defined  dependent 
variables,  i.e.,  variables  defined  at  a point. 
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4.5.2  Finite  Difference  Equations 


In  this  subsection,  the  numerical  scheme  for  computing 
an<l  updating  the  dependent  variables  defined  in  the  previous 
sections  is  described.  The  convention  used  here  defines  all 
quantities  at  integral  time  steps,  tn;  this  obscures  the 

centering  of  some  equations,  but  has  the  advantage  of  identi 
fying  time  t 


n 


with  cycle  number  n.  However,  the  equations 


re  not  entirely  time -centered ; such  centering  would  greatly 
increase  the  length  and  difficulty  of  numerical  procedures 
especially  the  pressure  iteration.  Since  there  is  only  one 
independent  component  of  the  Jeviatoric  stress  tensor  in  plane 
and  spherical  geometry  (using  the  appropriate  principal  axis 
coordinate  system) , but  two  in  cylindrical  geometry,  and  since 
cylindrical  geometry  is  a case  of  little  interest,  the  equa- 
tions have  been  written  for  only  plane  (d-1)  and  spherical 
(d=3)  geometry. 

The  following  is  a description  of  the  current  scheme 
used  for  setting  up  the  computational  grid,  and  advancing  the 
grid  through  one  computational  cycle.  No  attempt  is  made  to 
justify  the  centering  or  the  calling  sequence  of  the  equations 
hue  to  the  complexity  of  the  phenomena  treated,  the  only  adequate 
test  of  the  numerical  method  is  comparison  of  computations  with 
known  solutions.  Variables  will  be  defined  as  they  are  encoun- 
tered. The  space  centering  convention  used  is:  cell  boundaries 

and  the  tuff  and  water  velocities  are  defined  at  integral  space 
points  x.;  all  other  quantities  are  defined  at  cell  centers, 
1+1/2.  Note  that  subscripts  in  this  subsection  refer  to  spatial 

positions,  not  vector  and  tensor  components  as  in  the  previous 
subsection . 

In  the  lb TUP  routine  the  grid  boundary  array  x;  is  set 

together  with  tuff  mass  1 water  mss  tm)  . 

1 + 1/2’  Vclter  mass  + crushup 

parameter  a{l1/,,  and  water  volume  fraction  n?.,^.  The 


values  of 


1 + 1/2 
*1  i+1/2  and 


, ...  .i+1/2 

A2  i+1/2  are  lnitialized  to  unity. 
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The  updating  of  all  variables  in  a computational  cycle 
is  currently  done  in  the  following  sequence,  which  employs  a 
velocity  boundary  condition;  a scheme  for  incorporating  a 
stress  boundary  condition  is  under  development. 


First,  the  time  step  dt1 


criterion 


.n+1  _ . ( dxi+] 

- mm  jjy- 

( Ci?] 


' “ + a 
'i+1/2  7 


is  computed  from  the 


where 


i n 

dxi*l/2 


n n 

i+1  ' xi  * 


Nj+1  is  the  Courant  number  ( n£  + 1 < l)  , and  ^ and  V 
are  the  tuff  and  water  sound  speeds.  The  time  step  is 
the  minimum  value  of  the  indicated  quantities  over 
all  cells  in  the  active  grid  (i.e.,  all  cells  containing 
nonzero  stresses) . 

Next,  a test  is  made  to  determine  whether 
(1)n  „ , |d)n  | 

Pl-1/2  1 fa  * Mx  | Pi. 1/2  j (4.« 

where  I is  the  active  grid  counter  and  f is  an  accuracy 

factor  10  6);  if  this  condition  is  satisfied,  I is  advanced 
by  one . 

The  left-hand  boundary  of  the  grid  is  defined  to  be 
point  i=0.  Using  a velocity  boundary  condition,  one  specifies 


UV!  0)ntl 


where  f(t)  is  some  function  of  time.  The  left-hand  grid 


1 S 6 


boundary  is  then  updated: 


n + 1 n A ^n+1  ,.n+l 

*0  = x0  * v0  dt 


(4.67) 


It  is  assumed  that  this  updating  corresponds  to  a physical 
process  m^wJUch  worlds  done,  so  the  tuff  and  water  total 
energies,  W and  W , of  the  first  cell  are  updated 
according  to  the  prescription 


'J’n.l  . 

Vl/2 

CDn 

wih 

* £aI 

' n+1 

,xo 

r cvon+i 

dtn+1( 

pl/2 

q 1 / 2 

-T"  ) 

bl/2  / 

1 

'?■*!  . 
vl/2 

(2)n 
w 11 
" 1/  2 

* £a( 

>+1 

(Xo 

r 'v^1 

dtn+1( 

<«n 

pl/2 

. (2J*  ) 

ql/2/ 

(4 

.68) 

where  the 

area 

factor  is 

given  by 

fA 

= 1 

for  d = 

i . 

(4  . 

.69) 

fA 

* 4tt 

for  d = 

3 . 

CD 

(2) 

and  q 

and 

are 

the  tuff  and 

water 

artificial 

vis  - 

cosities,  defined  by 


<°>n 

qi+l/2 


(al 


for  d v 


(a) 

P 


n 


n Wn  / Wn  t“)n  \ 

1*1/2  d vi+l/2  (fqd  vi+ 1/ 2 " £<!  Ll*l/2 / 


i+1'2 


< 0 , 


(4.70) 


(°0n  (<On  («)n  (“)n  (°0n 

qi  + l/2  = "pi+l/2  d vi+l/2  fJL  Ci+l/2  for  d vi+l/2  - 0 * 


(1) 

S is  the  principal  deviatoric  stress  in  the  tuff. 
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I!ere»  fq  (-  1-6)  and  f£(^  0.2  5)  are  the  quadratic  and 
linear  artificial  viscosity  coefficients,  and 


(a) 


d v 


n 


i+i/: 


(a) 
v 'n 

vi  + l 


(a) 
v . 


n 


(4.71) 


Thus,  linear  and  quadratic  artificial  viscosity  terms  are 
used  when  a cell  is  crushing  up,  hut  only  the  linear  term 
is  retained  when  a cell  is  unloading. 


After  the  velocity  boundary  condition  is  applied, 
the  pressure  equilibration  conditions  are  applied: 


pi*l/2'  pi+ 1/2 1 


vn+l 


(2) 


i + 1/2’  ni  + l/2  ’ L i + 1 / 2 


n+1 


(1) 


n+1 


and 


^n+1 

1 + 1/2 


(4.72) 


are  computed  according  to  the  procedures  described  in  the 
previous  subsections.  Also,  in  the  pressure  equilibration 
routine,  the  tuff  and  water  energies  are  updated  by  setting 


(1) 

'Vi  + l/2  " 


(1) 


W. 


n 


i + 1/2 


+ AW 


n 

i + 1/2 


(4.73) 


(2)* 
lVi  + l/2 


(2)n 
w n 

"i+1/2 


AW1? 


i + 1/2 


where 


AW 


n 

i + 1/2 


/( 2) 

1 Pi 


nr 


n • 


n+1 


n+1 


i + 1/2 


(2) 

„ n 

qi+ 1/2 


W(2)n+1  (2)n  ) 

M ni  + 1/2  * ni  + 1 / 2 / 


i + 1/2 


(4.74) 
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where 


Vi  + l/2  = fv[(xi  + l ) ' (xi  ) ] 

and  the  volume  factor  is  defined  by 

fy  = 1 for  d = 1 , 

fy  = 4 m/ 3 for  d = 3 . 


(4.75) 


(4.70) 


Next,  intermediate  values  of  the  artificial  viscosities 
are  cajcylated  as  in  the  formula  above,  with  V”  replaced 
with  Ci+i/2-  Values  of  p 3i  are  calculated  by 


(D* 

p B-  = u 


/('n')n  + 1 + ('n')n+1 

5T  \ ni+l/ 2 ni -1/2 


f ('•!■  - "!") 


Xn  - xn 

1*1  xi-l 


, C2),  C2)ntl  (2), 

Pi-l/2  qi-l/2  * Pj.1/2  * qin-l/2 
" 1 1 — £ 


i-1/2 


TTJ 


) 


n . 


n+1 
i + 1/2 


/ (2) 

\ ni 


n+1 
i + 1/2 


^n+! 

ni  -1/2 


) 


(4.77) 


The  tuff  velocities  are  calculated  by 


(1) 

v. 

i 


n+1 


CDn 

V 

1 


+ m 


2 dt 


n+1 


m . 


i-1/2 


TO 

m i + 1 / 2 


H)  [pi- 


/ r 

ifA 


CD. 

qi+ 1/ 2 


n + 1 

1/2  + q i -1/2 


(1) 


n 


(1) 


n+1 


3 i -1/2  ‘ Pi  + 1/2 


+ + l k (CQ)n  j xi  + l ‘ xi-i] 

i + l/2  7 d \ i-l/2  Si  + l/2  / ^Ti J 


+ 7 £v[(xi  + l)d  - (xi-l)d]  p(B-*( 


(4.78) 
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I 


where 


5tl  = 0 for  d = 1 , 


= 1 for  d = 3 . 


(4.79) 


The  Tuler ian  coordinates  (cell  boundaries)  are  updated 
according  to 


x^1  - xj  * at" 


(4.80) 


(2) 

The  water  velocity  v is  updated  in  three  steps,  a 
procedure  necessary  for  momentum  conservation.  The  first 
step  is 


(2) 


<2>n 

vi  + JTJ 


2dt 


n + 1 


n 


m ■ - , , + 
l -1/2 


(2) 
n n 
Pi+l/2 


i + 1/2 


' f 

I A 


(*nc% 


1/2 


(2)* 
qi  - 1/2 


(2)*  \ 
qi  + l/2  / * 


1 f 

7 fV 


CD  j 

P Si  j 
(4.81) 

The  second  step  sweeps  over  interfaces  between  "momentum  cells" 
and  uses  the  so-called  "donor  cell"  method,  which  is  known  to 
yield  stable  results  in  hydrodynamic  calculations  involving 
advection,  i.e.,  transport  of  matter  between  cells.  At  each 
interface  between  momentum  cells,  corresponding  to  x.  , 
the  difference  between  water  and  tuff  velocities  at  time^t 
is  calculated:  n 


dvi+l/2 


- ^ . <J>n) 

7 \ i + I i vi  + l 'if  • 


(4.82) 
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Tf  dvi  + i/2  1 °*  mjmcntu^  ceH  i is  called  the  donor  cell; 
if  dvi  + i/2  < donor  cell  is  momentum  cell  i+1.  Then, 

at  each  interface  xi+i/2* 


(2)**  (2)*  2 /(2)(2U* 

vi  + l = Vi  + 1 + vrrv 1777“  6lm  V 

mi +1/2  + mi+3/2 

(2)  **  ( 2 ) * 2 /( 2)(2)V* 

v i = vi  ‘ tt; ctt“  6\ m v ' i+i/2 

mi -1/2  + mi  + l/2 


(4.83) 


whe  re 


* 

i + 1/2 


dt 


n*-l 


n (2>n  n 

- ‘5  V dv?.l/2 

where 

D = 

i or  i+1,  depending  on 

whether  momentum 

i or 

i+1 

is  the  donor  cell,  and 

1 /C2) n C2)n 

\ 

n 

n 

7 \ mn-l/2  + mD+ 1/2  / 

PD 

f 1 l"1(xn  + xn  \ld  - 1 

fv(l7\  ^ xn+i)J  [7 

(xd - 1 + xn)]  | 

The 

water  mass  in  each  cell 

is  updated  using 

donor 

cell 

method : 

(2)"*1  (2,n  , ,| 

(2)  \ * 

mi+l/2  ' mi+l/2  Sl 

m ) i > 

(2)"*i  C2)"  A 

(2)v* 

rai-l/2  * mi -1/2  • M 

m ). 

(4.84) 


(4.85) 


(4.86) 
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where 


■f™’)*  - at-1 


°n  Jv"  > 


(4.87) 


*5  • - “!■ , 


(4.88) 


n 

cn  = 


(2) 


m 


n 


:v  _(xD+l/2  ) * ( xn - 1/ 2 ) 


(4.89) 


Mere  n = i - 1/2  if  dv?  > Q,  and  D = i + 1/2  if 
dv?  <0. 

A final  momentum-conserving  correction  to  the  water 
velocity  is  made  using  th°  updated  water  masses: 


(2) 

Vi 


n+1 


(2) 

vr  • 
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** 


(2) 

„ n 

m . i / t + 

i-l/2 


(n) 


m . 


r?7 


_ n+l 
mi - 1/2 


TTJ 


i + l/2 


m 


n 


(4.90) 


i + l/2 


New  intermediate  values  of  the  artificial  viscosities 


(a)**  ,a, 

q i + 1 / 2 ’ are  calculated  using  the  above  formula  with  the  v^ 
replaced  by  v^  . Next,  the  cell  total  energies  are 
updated  by  sweeping  through  cell  interfaces  x^- 


( 1)  ** 
Wi+l/2 


(1)  ** 

W.  , 

l - 1/2 


CD* 

Wi+l/2  + 6 


(1)* 

^i - 1/2  * 6 


I (1) \ ** 

( * )i  - 

( " i 


(4.91) 
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-1 


where 


/ CD\ ** 

l W )i 


= A 


n+1  ,^n  + l CD 


dt 


vi 


n+1 


X 1 
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n+ 1 CD  ** 

1/2  + qi-l/2 


(1) 


S. 


n 


CD 


n+1 


i'1/2  + Pi+l/2 


CD** 

qi+l/2 


CDn 

S 11 
i+1/2 
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(4.92) 


For  the  water  energies, 
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(4.93) 


(4.94) 
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w 


'n  ■ fv 


(xn  \d  | 

l.n  \d‘ 

.1  D+l /2>  \ 

in- 1/2  ! 

(4.95) 


"ere  n = i-1/2  if  dvj*1  > 0,  n . in/2  if  dvj+1  < 0 . 

Next  in  the  computational  sequence,  the  A.  are 
updated  and  new  deviatoric  stresses  calculated. 


kn  + l 
1 i + 1/2 


n+i 

A2  i+1/2 
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(d=l) 


(4.96) 


,n+l 

A1  i+1/2 
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1 i + 1/2  I 1 + dt 
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n+i  vi+l  vj 
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CD*  CD  *t  m _ / n + 1 

l*l/2=  Si*l/2  (A1  iil/2  ’ A1  in/2  * A2  + i*l/2 
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Finally,  corrections  are  made  to  the  water  and  tuff  energies 
due  to  interaction  terms: 
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4 . 6 STRESS  PULSE  PROPAGATION 

The  thermodynamic  POROUS  code  has  been  exercised  in 
the  planar  option.  Initially,  a calculation  was  run  for  the 
following  input  parameters: 

0)  (2)  . , 

v = v = 5 x 10  cm/sec 
0 0 


pe  = 2.22  g/cc 
1 o 


pe  = 0.9982  g/?c 
2 0 


dx  = 0.02  cm 
0 


D = p/k  = 0.25  x 103 * * * * * 9/sec 


p =50  kbar 
P 


Yp  = 5 kbar 

p = 15  kbar 
1 c 

(1) 

n =0.8 


(2) 

n =0.1 


(3) 

n =0.1 

o 

The  results  of  this  calculation  arc  shown  in  Figs.  4.1  through 

4.3.  Figure  4.1  shows  the  particle  velocities  at  t = 3 psec. 

There  is  roughly  a 10%  overshoot  in  water  velocity  at  the 

shock  front.  This  result  is  qualitatively  similar  to  the 
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Fig.  4.2- -Volume  fraction  prcfils  for  the  tuff 
and  water  components  at  t = I nwr 


(I)yr:es /err  ■ 10’) 


1 


mechanical  POROUS  code  cul ations  reported  in  3SR-648. 

The  volume  fractions,  n and  n arc  plotted  in  Fig.  4.2. 

The  tuff  volume  fraction  experiences  an  increase  from  0.8 
to  approximately  0.908  through  the  shock  front.  The  water 
volume  fraction  changes  from  0.1  to  0.092.  An  interesting 
phenomenon  occurs  at  the  left  boundary.  Here  the  water 
volume  fraction  drops  from  0.092  to  roughly  0.083.  This 
results  from  the  advcction  of  water.  At  high  enough  stress 
levels,  this  can  become  a serious  problem  as  the  water 
volume  fraction  in  the  boundary  cells  may  approach  zero.  We 
arc  still  in  the  process  of  making  the  POROUS  code  general 
enough  to  handle  this  problem  as  well.  In  Fig.  4.3,  we  show 
the  axial  stress  profile,  o , at  t = 3 usee.  In  contra- 
distinction to  the  mechanical  POROUS  code,  no  overshoots  or 
oscillations  are  observed  at  the  shock  front.  This  is,  of 
course , the  result  of  using  a leap-frog  type  finite  difference 
scheme  instead  of  the  Lax-Wcndroff  scheme  employed  in  the 
mechanical  POPOUS  code. 

To  evaluate  the  effect,  of  porosity  and  water  content, 
two  additional  calculations  were  run  using  the  following 
initial  volume  fractions  (all  oti.er  parameters  were  kept  the 
same  as  in  the  above  calculation): 

(1) 

(1)  n = 0.9 

o 

(2) 

n =0.0 


(3) 

n =0.1 

o 

(1) 

(2)  n =0.9 

o 

(2) 

n = 0 . 1 


(3) 

n = 0.0 

o 
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The  stress  profiles  for  all  three  cases  at  t = 3 psec  arc 
plotted  in  Fig.  4.4.  The  dry  case  (1)  results  in  somewhat 
higher  Jjress  and^yve  velocj|y  than  the  partially  saturated 
case  x n^  - 0.8,  n^  = 0.1,  n^  = 0.1  j.  The  completely 
saturated  case  (2)  leads  to  still  higher  stress  and  wave 
velocity.  These  calculations  demonstrate  that  the  major 
effect  of  void  porosity  is  to  lower  the  stress  and  wave 
velocity  ampl i tu  les . 

These  calculations  are  of  a preliminary  nature.  The 
code  is  presently  being  modified  to  include  a pressure  boundary 
condition.  Also,  an  attempt  is  underway  to  generalize  the 
boundary  treatment  to  handle  high  stress  levels.  During  the 
near  future,  the  new  POROUS  code  will  be  exercised  in  spheri- 
cal geometry  as  well,  and  results  will  be  compared  with  those 
obtained  from  homogenized  models. 
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Fig.  4 . 4 - -Comparison  of  stress  profiles  at 
t = 3 usee  for  a porous  tuff  matrix  whose 
pore  volume  is  hal f - saturated , saturated, 
and  void  of  water. 
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ELASTIC  I LIJ I 0- ROCK  I NTF.RAO  I'  ION  AS  A MFiCllANISM 
i;OR  TRIGGERING  EARTHQUAKES 

5,1  illp.  I>0!  1 ^ AL  FOR  TRIGGERING  EARTHQUAKES  RY 

AHERIN'G  TllF,  GROUND  WATER  CONDI  T I QMS 

The  crust  of  the  earth  is  laden  with  tectonic  stresses 
that  heave  masses  of  earth  to  form  mountain  ranges  and  drive 
continents  to  move  one  with  respect  to  another.  Earthquakes 
tend  to  occur  in  those  regions  where  the  magnitudes  of  the 
principal  stresses  become  widely  different  (high  levels  of 
shear  stress).  This  very  intense  state  of  stress  is  generally 
associated  with  the  earthquake  belt  which  marks  the  boundary 
between  the  crustal  plates;  however,  as  evidenced  by  the  oc- 
currence of  past  earthquakes  throughout  the  world,  critically 
high  levels  of  shear  are  generated  over  a much  wider  region. 

In  many  regions  of  the  United  States,  tectonic  stresses 
are  approaching  an  unstable  condition.  This  condition  exists 
today  just  as  it  has  for  centuries.  In  the  recent  several 
years,  however,  the  possibility  of  triggering  an  earthquake 
by  artificial  means  has  become  a problem.  For  example, 
earthquake  activity  has,  in  some  cases,  increased  in  the 
vicinity  of  a newly  constructed  reservoir.  Following  the 
construction  of  Ko/na  Dam  of  South  India  and  the  subsequent 
filling  of  the  reservoir,  a sequence  of  small  earthquakes  were 
recorded  over  a period  of  several  years.  Then,  on  December  11, 

)7,  3 maSnitude  6.5  earthquake  occurred,  which  resulted  in 
considerable  loss  of  life  and  a major  economic  loss.f59] 

The  region  was  not  noted  for  earthquakes  prior  to  the  con- 
struction of  the  reservoir;  consequently,  it  appears  likely 
tiat  the  presence  of  the  reservoir  somehow  served  to  trigger 
the  catastrophic  release  of  stored  strain  energy. 

The  process  of  injecting  fluids  into  a tectonically 
stressed  region  is  another  means  by  which  earthquake  ruptures 
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can  be  artificially  initiated.  'lost  investigators  have  con- 
cluded that  the  swarm  of  hundreds  of  small  earthquakes  near 
Denver,  Colorado  in  the  period  following  1962  were  triggered 
by  the  pumping  of  waste  fluids  into  a 3671 -meter  disposal  well 
by  the  Rocky  Mountain  Arsenal.  A conceptual  model  was  pre- 
sented in  1968  by  Ilealy,  Ruby,  Griggs,  and  Rayleigh  which 
accounts  for  the  triggering  mechanism  of  the  injected  fluid. 

More  recently,  the  controlled  injection  of  water  "in  the  Rangley 
field,  Western  Colorado,  by  IJSGS  researcheis  ,62 ^ clearly  de- 
monstrated tne  potential  of  tuis  process  for  triggering  earth- 
quakes. In  this  experimental  program  water  was  injected  into 
the  western  portion  of  the  Rangely  Field  along  a previously 
mapped  fault  zone.  As  the  water  was  driven  into  the  faulted 
region,  the  fluid  pressure  mounted,  and  small  earthquakes  be- 
gan to  occur.  When  the  fluid  injection  was  stopped,  the 
occurrence  of  earthquakes  persisted  for  a period,  then,  as 
the  fluid  pressure  began  to  diminish,  the  earthquake  activity 
also  diminished. 

Actually,  this  description  of  the  earthquake  activity 
associated  with  the  fluid  injection  process  is  partially  based 
on  conj ecture , since  the  fluid  pressure  interspersed  in  the  rock 
was  not  actually  measured  at  the  point  cf  rupture.  The  sequence 
of  observed  activity  is  consistent  with  the  analytic  expressions 
developed  near  the  end  of  this  section  for  a simplified  spheri- 
cally symmetric  fluid  injection  system.  The  theory  indicates 
that  there  is  little  possibility  of  triggering  an  earthquake 
at  the  time  that  the  pumping  is  first  begun;  it  takes  some 
^ime  for  the  fluid  pressures  to  spread  over  the  subsurface 
region  of  incipient  rupture.  On  the  other  hand,  when  the 
pumping  is  stopped,  the  fluid  pressures  continue  to  grow  for 
some  time  before  leveling  off  and  diminishing  to  the  background 
level.  Hence,  based  on  our  theoretical  model,  we  would  expect 


i 
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the  earthquake  activity  to  continue  for  some  time  after  the 
fluid  injection  is  stopped.  The  time  lag,  of  course,  depends 
on  the  subsurface  porosity  and  the  regional  extent  of  the 
fault  zone. 

Archambeau  and  his  co-workers  at  the  California  Insti- 
tute of  Technology  are  just  beginning  a program  to  monitor 
earthquake  activity  associated  with  the  injection  of  fluid 
into  the  Santa  Fe  Springs  Field,  Los  Angeles,  We  anticipate 
a continuous  exchange  of  data  with  this  group  in  order  that 
we  might  test  our  theoretical  description  of  the  mechanical 
interactions  between  the  pore  fluid  and  the  rock  as  it 
relates  to  the  field  experiment. 
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5 . : vlh.rx  vn:  tih-orm  ical  loitHULAnoxs 

Various  theoretical  models  have  been  developed  to 

describe  the  mechanical  interactions  between  the  solid  and 

the  fluid  constituents  of  a saturated  porous  solid  material 

>uch  as  soil  or  rock.  Much  ol  this  work  has  been  mot  i tilted 

by  engineers  concerned  with  the  gradual  settlement  ol  satura 

ted  oils.  A simple  mechanism  to  explain  this  consolidation 

i i ...  . . [63] 

process  was  first  proposed  by  lerzigni. 

The  next  major  extension  to  the  theory  of  consolidation 
was  made  by  Biot[b4]  in  which  the  linear  consolidation  pro- 
cess was  modeled  in  three  spatial  dimensions.  Unfortunately, 

Biot's  formulation  does  not  explicitly  show  how  the  various 

volumetric  strains  enter  the  analysis.  He  states  explicitly 

that  the  pore  water  is  considered  to  be  incompressible.  We 

assume  from  this  that  the  small  grains  of  solid  material  would 

also  be  considered  to  be  incompressible.  On  the  other  hand, 

Biot's  theory  lias  proven  satisfactory  for  explaining  expeii 

mental  results  for  the  consolidation  of  a solid  sphere  of 

saturated  clav,t65,66]  and  for  the  consolidation  of  a two- 

r / 7 ] 

storied  aquifer. L ' J 

In  Biot's  formulation,  and  in  essentially  all  subsequent 
formulations,  we  find  two  physical  constants  that  serve  to  account 
for  interactions  between  the  two  constituents.  These  constants 
are  operationally  defined  through  explicit  laboratory  tests. 
Although  it  is  intuitively  obvious  that  *hese  terms  arise  from 
changing  dimensions  of  the  fluid-filled  pores,  no  effort  is 
made,  at  least  in  Biot's  formulation,  to  develop  the  physical 
processes  that  give  rise  to  the  two  terms.  Details  of  this 
type  are  needed  for  isolating  the  rock  stress  (as  opposed  to 
the  composite  stress)  upon  which  to  base  fracture  criteiia. 


5.3  QUASI  STATIC  LINI-AR  TINC  ASSUMPTIONS 


In  subsurface  geologic  formations,  where  stresses  are 
high,  the  compressibility  of  the  ground  water  and  the  rock 
grains  is  likely  to  have  a significant  effect  on  the  mechanical 
interactions.  lor  this  reason,  we  use  the  TINC  framework  to 
develop  equations  for  describing  the  mechanical  interactions 
between  porous  elastic  rock  and  ground  lluid  that  is  permeating 
through  the  pores  of  the  rock.  The  formulation  is  linear, 
and  as  such,  it  applies  to  saturated  soils  or  rock  where: 

1.  The  strains  are  small  compared  to  unity. 

2.  The  stresses  and  strain^  in  the  consti- 
tuents are  linearly  related  by  isotropic 
elastic  constants. 

3.  The  velocities  are  slowly  varying  so  that 
inertial  forces  can  be  neglected. 

4.  The  drag  forces  between  the  pore  fluid 
and  tiie  rock  matrix  are  linearly  related 
to  the  relat  ve  veloc  ty  between  the  pore 
fluid  and  the  rock  grains  [Darcy's  law). 

5.  The  rock  is  saturated.  For  the  case  where 
small  gas  bubbles  are  present  in  the  pore 
liquid,  we  assume  that  the  gar.  moves  with 
the  liquid  influencing  only  the  bu'l; 
modulus  and  density  of  the  gcs-liquid  mixture. 

For  the  case  where  the  roex  contains  small 
voids  that  are  not  evnnected  [isolated  from 
the  permeating  fluid),  we  use  the  term  pore 

to  refer  to  the  connected  pores  and  consider 
the  isolated  voids  to  be  homogenized  into 
the  rock  grains  thereby  influencing  only  the 
bulk  modulus  and  density  of  the  graii.s. 


3 . 4 CO.A S I. UV ATI  Q.\  liOUAT  I OX'S 


I lie  general  1I\'(  formulation  has  been  presented  in 
aSk':°"  aiui  3SR-  04  8 , as  uc  1 1 as  in  Section  I\  of  this  report. 
However,  in  this  section  we  restrict  our  attention  to  linear 
behavior,  and  as  such,  much  of  the  reduction  of  the  basic 
conseix at  ion  equations  to  their  final  form  for  computer  pro- 
cess inc  is  unique  to  this  section.  Consequently,  uc  have 
decided  to  briefly  redevelop  the  basic  TIXC  equations  in 
their  desired  form. 

The  mechanical  interactions  between  porous  geologic 
lock  and  interspersed  ground  fluid  arc  governed  by  the  con- 
sort at  ion  equations: 

Conservation  of  'lass 


Conservation  of  Momentum 


+ P t\  + p 8. 


Cons ? rvat i on  of  thermal  energy  can  he  disregarded  in  this 
uevc lopment ; also,  mass  transfer  due  to  chemical  interactions 
and  pnasc  changes  are  not  being  considered.  The  notation 
adapted  in  previous  TIXC  developments  is  being  used  with: 

(a)  = 1 for  the  solid  constituent 

= 2 for  the  fluid  constituent  fgas-liquid  mixture) 
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(*)e 
V c 


volume  of  constituent  (a)  excluding  the  volume 
occupied  by  the  complementary  constituent 

U)e  (2)e 

V 1 V (5. 

(a) 

V°A'  (5. 

volume  fraction,  hence 

(1)  (2) 

n + n = 1 (5. 

partial  density,  i.e.,  the  mass  of  consti- 
tuent (a)  per  unit  volume  of  composite 

composite  density 

CD  (2) 

P + 0 re 


velocity  of  a point  in  constituent  (a) 

(a) 

3 u . 

l 

~5F~ 


= displacement  of  a point  in  constituent 
(a)  from  its  starting  position  x 

partial  stress,  i.e.,  the  force  component 
in  constituent  (a)  per  unit  area  of  composite 

= composite  stress 

(1)  (2) 

= o . . + o . . rc 

lj  ij  U 

= body  force  per  unit  of  composite  mass  due 
to  gravitational  forces. 

= body  force  per  unit  of  composite  mass  that 
results  from  drag  forces  on  the  complemen- 
tary constituent. 

CD  (2) 

8i  * ■ 0 (5 
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5.5  CONST  1 run  Vi:  HQUAT  IONS 


In  order  to  complete  the  description  of  the  fluid-rock 


composite,  we  introduce  constitutive  equations: 

(1)  ( 2 ' 
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1.  I he  dray  forces  ( 
ve 


j = - ii  ^ are  related  to  tiie 
/elocitf  difference  between  the  fluid  and  the 
rock  v m through  an  extended  version  of 

Ua rev's  law 


(1) 


(2) 


D ^ = ■ 0 


C d 

o 


/(2)  (1)  V 

\ Vi  - Vi  ) 


(1) 

v 


(5.10) 


in  which 


d = - 


r~  r 

0 


(5.11) 


where  u is  the  fluid  viscosity  and  k is  the 
permeability  of  the  rock. 


2.  Changes  in  porosity  are  related  to  volumetric 
strains  in  the  constituents  by  a power  series 
expansion.  Retaining  only  the  linear  terms 
for  the  rock  constituent  we  write 


CD  (1)  (1)  (2) 

n / n = 1 + b c + b t 
0 1 2 


(5.12) 


(a) 


in  which  c is  the  volumetric  component  of 
partial  strain  defined  by 


Co) 

£ 


(o) 

* i j E ij 


(5.13) 


The  volume  fraction  n for  the  fluid  constituent 
is  obtained  from  Eqs . (5.5)  and  (5.12). 

3.  In  the  final  constitutive  equation  we  relate 
stress  and  strain.  First,  let  us  consider  the 
special  case  when  there  is  no  pressure  in  the 
pore  fluid.  Partial  stresses  in  the  rock  matrix 
are  related  to  the  partial  strains  in  the  rock 
matrix  by  the  conventional  isotropic  Hooke's 
law  of  linear  elasticity 


(1)  (1) (13  (1)  (1) 

Oij  = 2 y e ...  + A «...  e (5.14) 

(1)  (1) 

where  y and  A are  Lame's  constants.  An 
alternative  expression  is  obtained  by  ^compos  ing 
the  partial  sjjress  into  hydrostatic  oJ  and 
deviatoric  components 


(1)  (1)  (1) 

■ 0 hj  + sij 

such  that 


(1) 

o 


is.. 

3 ij 


CD 


a . . 

ij 


(1)  CD 

K e 


and 


(1) 

S.  . 
ij 


(1) 


CD 
«ij  0 


c1)  (ci)  . cm 

2y  \ e . . - ? 5. . e 1 
x ij  3 ij  / 


(5.15) 

(5.16) 

(5.17) 

(5.18) 

(5.19) 
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in  which  the  hulk  modulus 


(1) 

K = 


(1)  2 CD 

A + T * 


.Now  we  consider  the  case  in  which  the  pore  water  is 
acting  under  pressure,  h'c  will  assume  that  introduction  of 
poie  p i es  sure  does  not  affect  the  deviatoric  stress -strain 
lelatlonship  for  the  saturated  specimen,  l:q.  (5.19).  In 
order  to  develop  suitable  constitutive  equations  for  the 
hydrostatic  components,  we  introduce  the  concept  of  effec- 
ti\c  stiess  the  actual  stress  in  the  solid  rock  grains 
u\ oi aged  over  several  grains.  This  same  concept  of  average 
microscopic  stresses  also  applies  to  the  fluid  phase  so  that 
the  effective  hydrostatic  stress  is  simply  defined 


(°0C  ]_  (a) 

0 = TcTT  0 (5.20) 

n 

Ihe  corresponding  effective  volumetric  strain — the  microscopic 
volumetric  strain  in  a single  constituent  averaged  over 
several  pore  dimens  ions— is  kinematically  related  to  the  par- 
tial volumetric  strain 


(a) 

V 7 V 


(a) 

n 


(a)  / 

n 11  + 

o ' 


(5.21) 


consequently 


(a) 

e 


e 


(a) 


1 . 


(5.22) 


Linear  hydrostatic  stress -strain  behavior  for  the 
constituents  is^governed  by  the  bulk  modulus  of  the  isolated 
constituent  K e. 


(a) 


e 
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(5.23) 


or  by  applying  Eqs.  (5.20)  and  (5.22),  we  relate  the  hydro- 
static components  of  partial  stress  and  partial  strain  using 
a variable  volume  ratio 

(a)  (a)  (ot) 

a = n K e 

o 


(a) 

n 

(a) 
. n 


(i  * '?) 


(5.24) 


U) 

Substituting  n from  Eq.  (5.12)  and  retaining  only  the 
linear  terms,  we  get 


CD 

CD  (D e r 

CD  (2)1 

a = 

n K e (1  + 

b ) E + b £ ; 

Q L 

1 2 J 

and  similarly  from  Eq.  (5.5),  we 

get 

(2)  (2) 

CD  CD 

CD  CD  \ C 2) 

o = K e 

- n b e + 

1 - n - n b I e 

o i 

0 0 2 / 

(5.25) 


(5.26) 


We  invert  (5.25)  and  (5.26)  to  obtain  expressions  for  the 
hydrostatic  components  of  partial  strain 


and 


(1) 


CD  (2) 


a a 

nr  ' i r 


K 


(2) 


CD  (2) 

a a 

~rr  ir 


(5.27) 


(5.28) 
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whc  re 


r d)  i 

(1)  CD  CDe  [1  + >>  - nQ  d + bx  + b?) 

K = nQ  K f (TJ  1 ~ 

1 - n ,1  + b )J 

L o p J 


(K,e  r 

ii  = _L_  i + b 


:D  1 

n (l  + b + b ) 
o 1 2 J 


CD  1 

l + b - n (l  + b +b) 

1 0 1 2 J 


CDe 

R - -t- 1 ♦ b 

l + b 


(D  1 

n (l+b  + b ) 

n 1 9 J 


(5.29) 


U'e  havesuc^f)^ded  in  relating  the  bulk  modulus  of  the 
rock  matr  , K ( o = 9 ) , with  the  bulk  modulus  of  the  rock 
grains,  — e,  by  introducing  a mechanism  to  permit  variations 
in  the  volume  fraction,  Eq . (5,12).  Also,  we  have  arrived  at 
expressions  for  Biot's^64  ^ interaction  constants  H,  H ^ , and  R. 

In  this  linear  development  we  require  that  the  final 
state  be  independent  of  loading  path.  Therefore,  the  strain 
energy  density 


1 (DCD  1 (2)(2) 

4 e a + t e a 


(5.30) 


can  be  obtained  by  applying  the  load  in  two  stages: 

(2) 

1.  The  rock  matrix  is  loaded  with  a =0,  and  then 

CD 

2.  the  pore  pressure  is  applied  with  o = 0. 


The  combined  strain  energy  density  then  becomes 


1 (1)  (1)  (1)  (1)  1 (2)  (2) 

U = q-  £ a + e a + -*-  £ a 

i 1 1 2 1 i 2 2 


or,  if  we  reverse  the  order  of  the  loads,  we  get 

: fD  'D  , (2)  (2)  (2)  (2) 

U = £ t +tc  a + e o. 

ill  i 2 2 12 


where  the  subscripts  denote  the  load  stage. 

We  equate  the  two  expressions  for  strain  energy  density 
to  obtain  a particularized  version  of  Bette's  reciprocal  theorem 


(1)  (1)  (2)  (2) 


c a 

2 1 


£ a 
1 2 


(5.31) 


(2) 

a (1) 


using  Eqs . (5.27)  and  (5.28),  respectively, 
conclude 


From  this  we 


H = H 
i 


(5,32) 


and  consequently 


nr; 

K e 


b 

2 

tot; 


(5.32) 


Actually  b is  negative  and  b is  positive 
1 2 


1 


Thus  we  find  the  complete  stress-strain  relati onships , 
effective  and  partial  quant  it ies , ^ jyv^Jye  a total  of  five 

independent  physical  constants:  u , K (or  alternatively 

r c , ..  ,c  (f,)e  (.  1 J . , , 

K bv  Lqs.  (5.29)),  k , n , and  b (or  alternatively 

Q 1 

b by  Lq . (5.33)).  The  first  four  constants  are  non- interact i ve 
in  nature,  i.e.,  the  stresses  and  strains  in  the  rock  are  inde- 
pendent of  those  in  the  fluid  and  conversely  the  stresses  and 
strains  in  the  fluid  are  independent  of  those  in  the  rock  for 

the  special  case  when  b = 0 (consequently  b =0). 

1 2 

liquations  (5.25)  through  (5.28)  suggest  a number  of 
alternative  tests  for  obtaining  the  single  interactive  term 
b . For  example,  we  can  use  13q . (5.26)  to  deduce  b from 

1 (2)  i.(2)  v 

the  volume  of  water  e that  flows  from  a drained  (a  = 0) 
triaxial  compression  test 


:i)\  (i: 

n )/  n 

0 9 l 


(5.34) 


This  expression  leads  to  some  interesting  bounds  for  b . A 
one  extreme  the  freely  drained  specimen  is  compressed  but  no 
fluid  seeps  out  (k2)  c).  At  the  other  extreme  the 

volume  of  fluid  that  seeps  in  is  equa^yo  the  vqjume  that 
the  specimen  is  compressed  so  that  ( c = e,  e = -e). 

From  these  extreme  modes  of  behavior  we  arrive  at  the  bounds 


I U)  \ U)  / UJ  \ UJ 

- ( 1 - n )/  n (l  - n )/  n 

uj— or;  b i / uj-  Tm 

1 + K C/  K C ( 1 - K c/  K e) 


(D\  (1) 
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and  by  use  ot  Eq.  (5.33) , we  get 


U) 

/ (])\ 

- ii 

/ n 

1 - n 

' 0 / 

a 

r v\ — 

- b 

' o / 

(1) 
/ n 


/ K 


K °/ 


(5.36) 


In  general  we  might  expect  the  b (b  ) for  saturated  clay 
soils  to  lie  near  to  tlie  lower  (upper)  bound. 

A critical  test  of  the  suitability  of  the  non-dif fusive 
constitutive  equations  developed  above  has  been  provided 
through  a series  of  experimental  measurements  recently  re- 
ported by  Nur  and  ByerleeJ68^  They  derived  an  "effective 
stress"  law  for  eliminating  the  influence  of  pore  pressures 
from  their  test  data  of  the  form 


/ CD  CUeU2)e 

<0ij>NB  * °ij  ‘ l1  ' K / K C)  o e5lj  CS  .37) 

using  riNC  notation  where  <0ij>i\qj  i-s  their  effective  stress 
term  intended  to  correspond  to  the  equivalent  stress  in  a dry 
specimen.  They  presented  experimental  evidence  that  their 
effective  stress  law  accnunts  for  the  influence  of  the 
pressurized  pore  fluid  considerably  better  than  the  conven- 
tional effective  stress  law  which  simply  subtracts  the  fluid 
pressure  directly,  i.e.,  Eq.  (5.37)  with  4 e = «.  The  same 
data  was  processed  using  the  TINC  framework  in  which  the 
effective  stress  in  the  rock  is  taken  to  be  the  microscopic 
stress  in  the  rock  grains  averaged  over  several  grains.  The 
rather  small  bias  in  the  wet  data  that  remained  after  the  attempt 
to  remove  the  influence  of  the  pore  fluid  by  Eq . (5.37)  was 
essentially  eliminated  using  the  TINC  concept  of  effective 
stress.  This  rather  impressive  demonstration  of  the  use  of 
the  TINC  model  to  fit  stress -s train  data  is  presented  in  de- 
tail in  Appendix  I). 
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5 . b FLllIO  SFCPAGH 

V.'e  now  havre  the  basic  equations  with  which  to  describe 
the  elastic  behavior  of  the  fluid-rock  composite.  Using  the 
constitutive  laws  developed  above,  we  shall  proceed  to  ex- 
press the  conservation  equations,  Fqs.  (5.1)  and  (5.2),  in 
terms  of  the  fluid  pressure 

( “) e -1  (2) 

P = - ° = Th  a (5>38) 

(l  - n ) 

(1)  , (1)  d(V  \ 

and  components  of  displacement  u.  and  e = - — for 

l \ dx.  / 

tae  rock.  1 


First  we  will  investigate  the  conservation  of  mass 
equations.  From  the  kinematics  of  the  deforming  composite, 
we  equate  the  ratio  of  density  to  initial  density  with  the 
inverse  ratio  of  volume  to  initial  volume  and  write 


(a) 


(a)  \ - 1 

e ) 


Consequently 


and 


(a)  (a) 

1 B p Be 

r*7  5t  ■ at 

p 


(a) 

1 Bp 


P 

o 


(a) 
B e 
Bx  ■ 

l 


(5.39) 


(5.40) 


(5.41) 
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when  we  disregard  components  of 
performing  the  indicated  derivat 
substitutived  into  Eq.  (5.1)  to 
the  conservation  equation 


strain  compared  to  unity 
ives.  These  expressions 
give  the  alternative  form 


after 

are 

for 


(a)  (a)  (a) 

' ^ P j ^ C 

P 


(a) 

v. 


3 


(5.42) 


Since  3x . - 9t  ' Clearly  the  first  two  terms  cancel  to 
first  ordef  accuracy .from  which  we  are  led  to  conclude  that 

tll?<Jhlrd  t6rm  vi(  3 £ /3xi  ) is  higher  order  than  the  term 
3 e ^3t  ’ [The  fact  that  the  mass  conservation  equations 
are  automatically  satisfied  is  a consequence  of  the  manner 

in  which  we  related  partial  density  to  partial  strain  in 
Eq.  (5.39)] 


Momentum  is  conserved  quasistatically , 
eliminate  the  inertial  forces  from  Eq.  (5.2) 
equilibrium  expression  for  the  fluid  phase  (a 


i .e . , we 
to  obtain 
- 2). 


the 


(2) 

3 o 


+ 


(5.43) 


The  expression,  in  the  absence  of  inertial  forces,  remains 
time  varying  because  of  the  drag  forces  between  the  fluid 
and  the  rock,  which  have  been  related  to  the  relative  velocities 
between  the  two  constituents  by  Eq.  (5.10).  Equation  (5.43) 
is  reduced  to  a single  nonsubscripted  equation  by  differentiating 
with  respect  to  xA  and  summing  on  i 


.(2) 

V2  o + V 


(2) 
3 c 
3t 


(1) 
3 £ 
3t 


(5.44) 
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(2) 

strain  component  e 
by  Eq . (5.26),  and 
pressure  by  Eq.  (5.38). 


(2) 

is  expressed  in  terms  of  o and 
o is  in  term  related  to  the  fluid 


In  carrying  out  these  substitutions  we  reason  that 
l * J , fa) 

P 3 is  of  order  p /3t  by  Eq.  (5.12),  and 

P(3  e /St ) is  of  order  (p/  K e)(3p/3t)  by  hqs . (5.26) 
and  (5.38),  and  jp/  p e ) (?  p/3t)  is  of  order  ^e^(3p/3t| 

by  hqs.  (5.26)  and  (5.38);  consequently,  the  product 
d i f fe rent i at ion 


3 

3 1 


(1) 

n 


(1) 
3 n 
3t 


to  first  order,  reduces  to  the  form 


(5.45) 


Similarly  we  reduce 


P 


(1) 

3 n 


to  the  form 


3 

Sx . 

l 


(USp 
n 3x^ 


lollowing  these  substitutions,  Eq.  (5.44)  becomes 


(5.46) 


(5.47) 
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wiic  re 


(2) 


1 - 


(1)  1 

n (1  + b ) 

0 2-1 


The  elastic  compliance  of  the  rock  m^tyix  influences 
this  "seepage"  equation  through  the  term  3 e / 3 t , which 
originates  not  only  from  volumetric  fluxuations  in  the  poro- 
sity (ncnzero  b and  b'  ) but  also  from  the  fact  that  the 

1 2 

drag  forces  depend  on  the  relative  velocity  between  the  rock 
and  the  fluid  rather  than  the  fluid  velocity  alone,  Hq.  (5.10). 
It  was  discovered  through  actual  computer  applications  that 
presence  of  the  interactive  term  3^e  /3t  strongly  influences 
the  seepage  process  to  the  point  where  rather  severe  stability 
problems  arise  in  some  cases.  This  troublesome  phenomenon  can 
be  eliminated  if  we  rearrange  Eq.  (5.47)  to^obtain  a modified 
interaction  term  as  some  combination  of  3 e /3t  and  3p/3t. 
The  modified  interaction  term  should  have  little  influence  on 
the  seepage  process  for  optimum  computer  processing.  The 
total  hydrostatic  stress  a * has  the  desired 

character.  As  ground  water  is  driven  through  the  subsurface 
rock  formations,  increases  in  the  pore  fluid  pressure  re- 
duce compressive  loads  on  the  rock  matrix  so  that  the  total 
hydrostatic  stress  tends  to  remain  invariant.  In  fact,  we 
find  that  the  total  hydrostatic  stress  remains  totally  in- 
variant when  fluid  is  injected  into  a spherically  symmetric 
environment,  Section  5.9. 


(1) 

a 


Sub^jtuting  from  Eqs.  (5.27),  and  (5.38)  with 
o - o , we  get 
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* i1  • + ft) 


(5.48) 
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waich  ue  use  to  eliminate 
arrive  at 


(1) 

3 c /3t  from  l:q . (5.4  7)  and 
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(5.49) 


1 


5."  ELASTI^J^^matjq^  in  THE  rock  matpty 


Equilibrium  conditions  in  the  elastic  rock  matrix  are 
expressed  most  conveniently  by  combining  the  momentum  equa- 
ons  or  tie  fluid  and  the  rock,  Eq . (5.2)  with  a - 1 and 

again  d y/llminating  ^ ^ f°rCeS’  E*-  (5.10).  We 
gain  disregard  inertial  forces  and  write 


3a.. 

n 


C2) 
3 a 


37“  + Pfi  = 0 . 
1 


(5.52) 


in  which  £.  = ^ 

f5  l 3 CD  i ' 'fe  USe  E<<s-  fS.151  (5.19),  Wd 

(o.27)  to  reduce  a into  counts  of  J'C f 

The  partial  stress  in  the  fluid,  is  then  ^ d * ’ 

pore  pressuie  using  Ea  f5  381  ana 

for  the  composite  becomes  equilibrium  equation 

2 4-  ((;us1). ) ♦ _i_  f(f«  2 onaf 

9 j ' 8xi  L K ' 7 u ' e J + pfi 


or  by  Eq.  (5.13) 

r,../  (n  to 


•:>  k [(■  • 


C 1 ) / 3 u.  3 U.  \ f /m  v Cl)  -] 

- y l-g~  + 2 (1)\  3 a,; 

j i \ ^T'J  ^7  y / -sxtj + 

/ d)\  , f/  Cl)  V 1 
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(5.53) 


(5.54) 
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For  the  special  case  of  homogeneous  material,  13q.  (5.54) 
assumes  the  foi  . 


or  alternatively 


(5.55) 


(5.56) 
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5.8  PR  I com. 


The  mechanical  interaction  between  ground  fluids  and 
poious  subsurface  formations  have  been  mathematical  1/  modeled 
using  a linearized  version  of  TINC.  The  penetration  of 
elastically  compressible  fluid  through  the  pores  of  a sub- 
surface geologic  mass  is  modeled  by  Eq.  (5.47)  or,  in  a form 
more  suited  for  computer  processing,  by  Eq.  (5.49).  The 
elastic  compliance  of  the  porous  material  to  the  pressure  of 
interspersed  fluid  is  represented  by  the  presence  of  an  in- 
teractive term  related  to  the  volumetric  strain  rate  in  the 
rock  matrix.  The  interactive  seepage  equation  has  the  form 
of  a diffusion  equation,  which  can  be  treated  by  one  of  several 
existing  numerical  codes.  A 2-0  finite  element  code,  ori- 
ginally programmed  by  i.'ilson , [ 69]  was  selected  to  treat  the 
seepage  process  with  the  term  that  contains  the  elastic  com- 
pliance of  rock  matrix  appearing  as  a source  term. 

I he  elastic  deformations  in  the  rock  matrix,  expressed 
by  Eq.  (5.54),  were  also  treated  by  a 2-D  finite  element  code 
in  which  the  pore  pressure  enters  as  a body  force.  The  two 
finite  element  codes  were  merged  into  a single  2-D  code  (FRI) . 
Figure  5.1  illustrates  the  sequence  of  operations  that  are  per- 
formed to  simulate  the  interactive  seepage  process.  The 
diffusion  equation  is  solved  to  give  the  fluid  pressure  field 
at  an  advanced  time  step  with  the  interactive  term  extrapo- 
lated from  the  previous  time  step.  The  updated  fluid  pressure 
field  then  feeds  into  the  elasticity  portion  of  the  FRI  code 
to  generate  an  updated  displacement  field  from  which  rock 
dilatation  and  total  hydrostatic  stress  are  computed.  The 
solution  to  the  diffusion  equation  is  repeated  to  give  a 
corrected  pressure  field  at  the  advanced  time  step  base  on 
an  updated  interactive  term.  The  complete  cycle  is  repeated 
two  to  five  times  for  a single  time  step  in  order  to  avoid 
any  lag  in  the  interactions  between  the  fluid  and  the  rock. 
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1'he  newly  formed  PRI  code  was  first  tested  on  non- 
lnteract ive  systems.  The  rock  elasticity  mode  was  suppressed 
o produce  purely  diffusive  seepage  which  was  compared  with 
analytic  solutions.  The  converse  situation,  in  which  the 
Iffusion  mode  was  suppressed,  was  also  run  to  test  elastic 
e formations  in  the  rock  matrix  against  closed  form  results. 

The  iist  interactive  test  was  performed  on  a 2-i) 
simulation  of  a 1-1)  compaction  process.  A closed  form  solu- 
tion was  oatained  for  the  response  of  a single  laver  of 
fluid  saturated  material  to  a uniform  surface  load.  Pig  5 2 
Plastic  deformations  in  both  the  rock  and  the  fluid  are 

treated  using  the  linear  version  of  TIMC  to  yield  the  series 
solution 


p(y)  = 


4p , 


i = 0 


(-1)) 
(21  + 1) 


e 


(2i  + l) 2 


it2  r/4 


sin 


(2i+l) 


(5.57) 


Pd  yo[HT  * TT-nHWu)] 

(dimensionless  time) 

d = fluid  diffusion  coefficient 
(large  d makes  seepage  slow) 

P - mass  density  of  the  saturated  composite 

n = porosity  of  the  solid  material  (the  ratio  of 
the  pore  volume  to  the  total  volume) 

A,u  = Lame's  constants  for  the  solid  material 

“ = bulk  modulus  of  the  fluid 


Finite  Element  Configuration 


Exact  series  solution 

• • Finite  element  solution 
t Dimensionless  time 

p Initial  fluid  pressure 

o 

y depth  of  impervious 
0 rock  layer 


Fig.  5. 2 --Comparison  of  FRI  code  solution  with  exact 
solution  for  a test  fluid/rock  interaction  problem. 


0 


5 = n 7 '(l'-'nTQV^T 

w = uniform  surface  load  density 


As  seen  in  fig.  5.2,  the  numerical  finite  element  solution 
rollows  the  exact  series  solution  well  even  up  to  late  times 
when  much  of  the  seepage  lias  taken  place.  A variable  time 
step  was  used  to  achieve  the  finite  element  results  with 
At  = d . 0 0 2 5 initially  and  At  gradually  increasing  to 
At  = 0.05  for  the  late  time  calculation  t >_  0.10. 

The  elastic  compliance  of  the  solid  material  as  the 
fluid  seeps  to  the  surface  follows  the  equation 


(1) 

£ 


np 

rr"nfr~^T 


(5.58) 


where 

(1) 

e = dilatation  of  the  elastic  rock  matrix, 

C1) 

= (1-n)  (xV2u)' 

(final  dilatation  state  when  p = 0) 

The  finite  element  calculations  follow  this  relationship  pre- 
cisely for  the  initial  est  problem. 
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5.1)  INJECTION  WI-LL 


'lore  recently  the  FRI  code  has  been  exercised  to  examine 
tlie  interactive  seepage  process  associated  with  an  injection 
uell.  Simultaneously  an  effort  has  been  made  to  obtain  a 
closed  form  solution  for  the  injection  of  fluid  into  a homo- 
geneous formation  with  which  to  compare  computer  generated 
results.  Obviously  an  analytic  expression  describing  this 
process  would  have  lar  reaching  implications  for  providing 
guidance  in  future  injection  projects.  This  effort  has  met 
with  considerable  success. 

Ue  consider  the  problem  of  a cased  well  (which  we  will 
ignore  in  our  analysis)  pumping  fluid  into  a spherical  cavity 
ot  radius  r^  which  is  located  in  a homogeneous  saturated 
region,  Fig.  5.3.  For  simplicity  we  assume  no  initial  stress 
in  the  rock  and  no  initial  pressure  in  the  fluid.  The  injec- 
tion process  is  then  accomplished  by  stepping  the  pressure 

in  the  cavity  instantly  from  zero  to  p where  it  remains 

o 

throughout  the  injection  process.  Pressures  develop  in  the 
region  outside  the  cavity  at  the  instant  the  cavity  pressure 
is  applied  (t  = 0).  These  stresses  can  be  calculated  from 
hqs.  (5.47)  and  (5.54)  with  d - 00  to  prevent  seepage  at 
this  initial  loading.  From  Fq.  (5.47)  we  get 


1 - n 


0 
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Cased  IVell 


I-ig.  5 . 3- - Spheri  cal  ly  symmetric  fluid  injection 
system.  The  pressure  is  introduced  in  the 

cavity  at  the  initial  time  and  then  held  at  this 
level  with  the  pore  fluid  in  the  rock  formation 
assigned  zero  initial  pressure. 


(1) 

where  K 


and  !1  are  defined  by  Eq.  (5.29)  and  Q is 
given  following  Eq.  (5.47). 

U'e  substitute  this  expression  for  the  initial  pressure 
p into  the  elasticity  equation,  Eq . (5.56)  to  get 


(1) 

u 
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)■ 


(5.60) 


using  the  fact  that 


(1)  CD 

V • U = £ 

> 


(1) 

Vx  u =0  (5.61) 


for  the  homogeneous  spherically  symmetric  injection  environ- 
ment where  only  radial  displacements  occur.  The  benavior  of 
the  rock  is  constrained  to  satisfy  the  boundary  condition  at 
the  well  injection  cavity 

c = - p at  r = r (5.62) 

rr  *o  o v ’ 

i ^ 

i • W * ) 

(1)  (1) 

c = - n p at  r = r (5.65) 

rr  o o o 


and  the  condition  of  no  displacement  at  large  distances  from 
tiie  well.  The  solution  to  Eq.  (5.60)  that  satisfies  these 
conditions  is  found  to  be 


(1) 

Uf  (x, 
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(5.64) 


which  leads  to 


(1) 

°ij (x,  t » 0) 


CD 
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(5.65) 


Somewhat  surprisingly,  we  find  that  the  initial  pressure  in 

the  cavity  generate,  no  hydrostatic  stress  in  the  surrounding 
rock  formation,  j e., 


U)  _ 1 (1) 

° 3”  i j °i j ~ ^ (5 .66) 


and  consequently  there  is  no  pressure  generated  in  the  pore 

fluid  outside  the  cavity  at  the  initiation  of  the  injection 
process. 

The  results  presented  above  apply  in  the  initial 
injection  process,  before  fluid  seepage  begins.  Throughout 
the  li °f  injection  well,  the  composite  stress 
°ij  °ij  + aij  must  satisfy  the  equilibrium  equation 

9a . . 

LL  , n 

0 (5.67) 

outside  of  the  injection  well  and  Eq.  (5.62)  at  the  cavity 
boundary  which  yields 


and  consequently 


3x . x . 
— 1 J 
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0 


(5.68) 


(5.69) 
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for  all  time.  This  simply  indicates  that  the  h>drostatic 
stress  in  the  combined  fluid- rock  composite  permits  us  to 
decouple  the  seepage  equation,  Hq.  (5.49),  from  the  elasti- 
city equations  for  the  rock,  since  -r-2-  --  0 for  all  time. 

Ihus,  our  "interactive"  seepage  process  is  described  by 

C V2  p = (5.70) 


u here 
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1 

p»d  (r  " nr  * n) 
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(1) 


(5.71) 


The  constitutive  terms  R,  K , and  I!  are  given  by  Hq.  (5.29). 

\t  the  initial  time  there  is  no  fluid  pressure  outside  the 

injection  cavity  as  deduced  above.  At  the  cavity  boundary 

the  pressure  is  held  at  p . From  these  conditions  we  are 

a 

able  to  uniquely  express  the  complete  fluid  pressure  time 
iiis  tory 


P(r,t) 


- erf 


/(r/V 


2/7 
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(5.72) 


where  the  dimensionless  time  t is  given  by 


o 


(5.73) 


Ke  substitute  this  result  back  into  the  elasticity  equations 
for  the  rock  to  determine  the  stresses  and  deformations  that 
are  occurring  there,  and  we  find 
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We  have  not  carried  out  the  indicated  integration  to  deter- 
mine the  radial  displacement  field;  however,  we  have  obtained 
the  early  and  the  late  time  limiting  values 
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It  appears  that  there  is  much  to  be  learned  from 
examining  the  results  developed  above  about  the  potential  for 
triggering  an  earthquake  by  pumping  fluid  into  the  eart'.  We 
note,  for  example,  that  shear  stresses  are  due  entirely  to 
the  presence  of  pressure  in  the  cavity;  there  are  no  shear 
stresses  generated  by  fluid  seeping  through  the  rock.  The 
greatest  potential  for  triggering  an  earthquake,  therefore, 
iies  in  the  tensile  stresses  (relative  to  the  pre -inj ection 
state)  that  are  generated  in  the  rock  matrix  as  pore  fluid  is 
pressurized  at  points  away  from  the  cavity.  Figure  5.4 
illustrates  the  spread  of  the  fluid  pressure  through  the 
pores  of  the  rock.  Both  the  time  rate  of  loading  in  the  rock 
matrix  and  the  pore  fluid  are  depicted. 

The  spherically  symmetric  injection  system  was  also 
treated  using  the  numerical  Fill  code.  Fluid  pressure,  rock 
stress,  and  rock  deformation  are  generated  with  the  FRI  code 
operating  in  the  completely  interactive  mode,  i.e.,  with  no 
assumptions  as  to  the  nature  of  the  interactions.  The  numeri- 
cal simulation  results  in  hydrostatic  stresses  in  the  fluid- 
rock  composite  which  are  two  orders  of  magnitude  less  than 
the  radial  and  hoop  stresses  in  the  rock  matrix.  Ideally  the 
hydrostatic  composite  stress  should  be  zero  for  the  spherically 
symmetric  environment,  Eq.  (5.69).  The  fluid  pressures  from 
the  numerical  simulation  are  presented  in  Fig.  5.4.  We  see  very 
good  agreement  with  the  analytic  solution  even  at  early  times 
where  the  pressure  front  is  quite  steep. 
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We  will  now  investigate  what  happens  when  the  injec- 
tion process  is  terminated.  At  time  t ^ (dimensionless 
time  Tj)»  the  pressure  in  the  injection  cavity  is  dropped 
l'rom  p to  zero.  This  can  be  accomplished  analytically  by 
simply  introducing  an  additional  injection  solution  with  a 

negative  cavitv  pressure  -p  at  the  retarded  time  t . 

o 1 

Then,  after  pumping  is  stopped,  we  get 


from  Eqs . (5.68),  ( 5.72),  (5.74),  and  (5.75),  respectively. 

We  find  from  Eq.  (5.80)  that  the  fluid  pressure  will  continue 
to  increase  away  from  the  injection  cavity  for  some  time 
after  the  injection  pressure  has  dropped  to  zero.  Hence  the 
potential  for  triggering  an  earthquake  is  not  eliminated  when 
the  pumpir^is  stopped.  Actually  the  most  critical  conditions 
(maximum  o^)  will  generally  occur  after  the  pumping  has 
stopped . 
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5-10  SUMMARY  AND  CONCMJST DM.q 


The  processes  by  which  alterations  in  the  ground  water 
state  can  trigger  a major  rupture  in  the  earth  should  be  care- 
fully examined  in  order  to  guard  against  the  possibility  of 
accidently  initiating  a catastrophic  earthquake.  The  possi- 
bility of  controlling  the  natural  process  should  also  be  con- 
sidered. It  might  be  that  earthquakes  can  be  produced  for 
the  purpose  of  relieving  stress  in  the  earth  that  might 
otherwise  accumulate  and  result  in  major  damage. 


The  theoretical  formulation  of  the  mechanical  inter- 
actions between  the  fluid  and  the  rock  is  somewhat  involved. 
Actually,  the  resulting  equations  are  quite  similar  to  those 
developed  previously  by  BiotJ64!  The  major  points  where  the 
linearized  TINC  formulation  differs  from  Biot's,  and  essen- 
tially all  subsequent,  formulations  are  noted  below. 


1.  A mechanism  for  fluid-rock  interaction  is 
provided  in  the  linearized  TINC  formulation 
by  dealing  with  variable  pore  dimensions. 
Using  this  approach  the  interaction  terms 
take  on  a new  meaning. 


2.  The  linearized  TINC  formulation  includes  the 
elastic  deformations  in  the  rock  particles  and 
in  the  pore  fluid. 


3.  The  rock  velocity  as  well  as  the  fluid  velocity 
is  considered  when  expressing  the  drag  forces 
between  the  fluid  and  the  rock  (Darcy's  law) . 

A finite  element  code  (FRI ) has  been  developed  to 
treat  the  TINC  seepage  equations  in  two  spatial  dimensions 
(plane  strain  or  axisymmetri c) . While  the  FRI  code  is 
suited  for  complex  geologic  formations,  we  have  concentrated 
on  elementary  cases  in  order  to  examine  the  accuracy  of  the 
numerical  simulation.  We  find  that  the  FRI  code  is  well 


209 


suited  for  interaction  calculations,  at  least  for  the 
cases  that  were  run.  Good  accuracy  is  obtained  even  one 
time  step  into  the  calculations,  Fig.  5.4.  The  FRI  code 
has  not  been  critically  tested  for  simulating  late  time 
solutions . 

Analytic  solutions  were  obtained  for  an  injection  well 
in  a spherically  symmetric  environment.  The  injection  well 
analysis  points  out  the  potential  for  triggering  an  earth- 
quake by  pumping  fluid  into  the  ground.  Whereas  the  stresses 
that  result  from  an  artificially  applied  load  (e.g.,  a 
surface  load  or  a pressurized  cavity)  diminish  as  1/r3 
away  from  the  point  of  application,  the  stresses  that  are 
generated  in  the  rock  matrix  by  fluid  injection  diminish  as 
1/r  at  late  times  (steady  state).  Furthermore,  the  rock 
stresses  generated  by  fluid  injection  are  not  relieved  by 
stopping  the  pumping,  in  fact,  the  stresses  continue  to 
mount  for  a period  of  time  after  the  pumping  has  stopped 
before  the  fluid  pressure  begins  to  diminish.  This  phenomena 
should  be  carefully  considered  in  order  to  avoid  triggering 
a major  earthquake  unintentionally. 
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VI.  DISCUSSION 


The  successful  implementation  of  the  TAMEOS  equation  of 
state  for  homogenized  rock-water-void  mixes  as  a subroutine 
in  the  SKIPPER  code  is  an  important  milestone.  This  subroutine 
can  be  readily  incorporated  into  other  standard  ground  motion 
finite  difference  codes,  whether  ID,  2D,  or  3D.  The  table 
look-up  routine  used  with  the  PEQ  model  can  also  be  used  in 
conjunction  with  tables  generated  with  the  P*EQ  and  PTEQ  models. 
This  capability  is  currently  being  added  to  TAMEOS. 

At  present  the  PEQ,  P*EQ  and  PTEQ  models  are  limited  to 
hydrodynamic  pressures  less  than  1 mbar,  with  primary  emphasis 
on  pressures  less  than  200  kbars.  Additional  work  is  needed 
to  extend  the  range  of  the  tables  to  be  used  in  TAMEOS  up  to 
tens  of  megabars.  Shock  metamorphisms ^ ^ in  the  rock 
component  (e.g.,  poreless  tuff  or  granite)  should  be  consi- 
dered to  determine  adequate  treatments  of  these  phase  changes. 
Previous  work  in  rock-water  mixes  by  Butkovich  ^ treated 
only  the  high  pressure  states  for  saturated  rock  under  the 
PEQ  model.  Cratering  calculations  are  believed  to  be  very 
sensitive  to  the  model  used  for  energy  partition  between  the 
components  upon  release  after  shock  processing. 

The  irreversible  pore  collapse  model  used  in  TAMEOS 
is  based  on  the  disconnected  pores  postulate.  Because  of 
the  sensitivity  of  the  stress  wave  propagation  calculations 
to  the  crushup  constitutive  model  at  these  low  pressures,  a 
study  should  be  made  to  determine  the  conditions  under  which 
the  postulate  is  appropriate.  The  connected  pore  model  and 
such  phenomena  as  partial  void  recovery  upon  release  should 
be  considered  in  conjunction  with  the  available  experimental 
data . 


The  inclusion  of  four  generalised  plasticity  models 
into  SKIPPER  provide  the  code  with  options  that  possess  both 
sophistication  and  considerable  flexibility  to  match  available 
laboratory  test  data.  The  cap  model  and  three  Mohr-Coulomb 
models  are  included,  one  without  work  hardening,  one  with 
isotropic  work  hardening  and  one  with  kinematic  work  hardening. 
Comparison  of  calculations  with  field  measurements  in  granite 
show  the  'lohr-Coulomb  model  with  kinematic  work  hardening  to 
give  better  agreement  than  is  possible  with  the  cap  model. 

It  appears  at  this  time  that  the  four  models  in  the  SKIPPER 
code  are  a sufficient  base  for  treatments  in  which  the  medium 
is  considered  as  a single  material.  Modifications  are  re- 
quired, however,  to  adequately  account  for  pore  pressure 
effects  and  relative  motion  between  blocks  of  jointed  rock 
masses.  These  phenomena  are  believed  to  be  the  basic  reason 
for  the  reduction  in  the  laboratory  flow  stress  value  (by  a 
factor  of  six)  that  is  required  to  bring  calculations  for 
granite  into  agreement  with  field  data.  Such  discrepancies 
do  not  appear  to  be  severe  for  competent  sedimentary 
materials . 

In  formulating  the  generalized  plasticity  models  for 
the  SKIPPER  code  a logrithmic  definition  of  strain  was  employed. 
This  is  a preferred  definition  for  one -dimensional  codes  but 
its  interpretation  in  2D  codes  is  not  apparent.  V.'ork  is  needed 
to  permit  these  plasticity  models  to  be  also  used  in  other  codes, 
e.g.,  CRAM. 

The  TINC  model  for  treating  relative  motion  between  the 
rock  and  water  components  in  ground  shock  calculations  has  been 
significantly  advanced  with  the  development  of  the  new  PORO'JS 
code.  The  ground  motion  predictions  for  a partially  saturated 
wet  tuff  using  this  new  thermodynamic  version  of  the  code  should 
be  modified  to  treat  a spherical  high  explosive  charge  as  energy 
source  so  that  the  comparison  can  be  made  with  an  instrumented 


L. 


212 


high  explosive  field  test.  Consideration  should  also  be 
given  to  modifying  the  TINC  plasticity  treatment  in  POROUS 
to  include  a generalized  cap  model  for  the  rock  component. 

The  2D  quasistatic  FRI  code  represents  a new  tool 
for  evaluating  the  interaction  between  a pore  fluid  and  a 
stressed  rock  matrix  as  the  fluid  is  driven  through  the 
rock  mass.  Ihe  code  should  be  applied  in  a series  of  cal- 
culations for  a region  around  an  injection  well  in  an  effort 
to  explain  the  mechanisms  associated  with  hydro  fracture . 

Ihe  intent  would  be  to  perform  the  quasistatic  analysis  at 
various  stages  of  rupture  in  an  effort  to  follow  the  redis- 
tribution of  tectonic  stresses  and  alterations  in  the  per- 
meability as  the  rock  is  fractured.  Calculations  should  be 
devised  for  comparison  with  field  data  on  the  Rangley  field, 
Colorado,  and/or  the  Santa  Fe  Springs  field,  Los  Angeles. 

Consideration  should  be  given  to  the  development  of 
a 2 - D dynamic  FRI  computer  code  with  a mechanism  for  spon- 
taneous rupture.  Such  a code  would  provide  a more  adequate 
treatment  than  currently  exists  for  evaluating  the  possibility 
of  triggering  earthquakes  at  NTS  by  the  passage  of  a shock 
wave  over  a pre-existing  fault  zone. 
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APPENDIX  A 

SOUND  SPEED  IN  A POROUS  MATERIAL 

The  model  of  a porous  material  introduced  in  3SR-648 
specifies  that  the  bulk  hydrodynamic  pressure  in  the  porous 
material  is  given  by 


where 


V = bulk  specific  volume 
E - specific  internal  energy 
a = distension  ratio 


(S3  (V) 

- n + n 

= m — 


n 


(S) 

n = volume  fraction  of  solid 


(V) 

n = volume  fraction  of  voids 

-(V  * E)  is  the  hydrodynamic  equation  of  state  of  the  solid 
material  where  Ve  is  the  effective  specific  volume  of  the 
poreless  solid,  i .e . , 


Ve 


(A. 2) 


In  the  present  formulation,  we  recognize  that  the  bulk  pressure 
is  averaged  over  the  total  surface  area,  so  that  the  effective 
stress  is  reduced  by  a factor  of  l/a.  This  differs  from 
Herrmann's  model  wherein  bulk  pressure  is  set  to 


^Herrmann  " -(a  ’ E)  ” £(Ve,  E) 


(A. 3) 


Preceding  page  blank 
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Carroll  and  llolt^^  have  discussed  these  twc  'ormulations 
and  recommended  that  l:q  . (A.l)  be  adopted  for  strcngthless 
porous  matrices.  They  did  not  rederiiC  the  expression  for  the 
sound  speed  of  the  bulk  material  in  their  analysis.  This 
relationship  is  presented  below  and  compared  to  the  original 
expression  proposed  by  Herrmann. 


Let  us  retain  Herrmann's  form  for  the  distension  ratio, 


1 . e . , 

a = a (p) 

Hie  sound  speed,  c,  may  be  obtained  as  10 ! lows 


£l  = -/i£\ 
V2  I'Vip 


+ P 


Oh  }y 


(A.  4) 


(A.  5) 


l)i  f ferenti  ating  lq.  (A.l)  leads  to 


and 


(A. 6) 


(A.  7) 


Since  tnc  sound  speed  in  the  solid  material  is  defined  by 


c 


2 

s 


3P  BP 

+ V2  P r 

3VC  C - 3P  ’ 


(A.  8) 
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we  arrive  at  the  following  expression  for  bulk  sound  speed 


i + (ve  -z_  + P \ allEl 

V 9Ve  -)  a 2 

In  the  limit  of  P -*■  0 , this  reduces  to 


(A. 9) 


(A. 10) 


where 

K = normal  bulk  modulus 


CS o = sound  sPeed  in  solid  material  at  zero  pressure 


Herrmann's  expression  is 


(c2) 


Herrmann 


(A. 11) 


It  is  apparent  that  both  results  indicate  that  major  reduc- 
tions in  the  sound  speed  of  porous  materials  can  be  achieved 
only  if  the  matrix  deforms  under  infinitesimal  compressions 
(ot  (0)  < 0).  i tie  major  difference  between  the  two  models  is 
the  factor  of  which  appears  in  the  numerator  of  Herrmann's 

expression.  This  term  is  troublesome  from  a physical  view- 
point since  it  implies  that  in  the  limit  of  a perfectly  rigid 
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mat  rice  structure  (a'(0)  -<•  0),  the  sound  speed  is  increased 
by  adding  to  the  void  volume  fraction.  In  the  present  for- 
mulation, the  a factor  no  longer  appears.  Hence  in  the 

o 

same  limit  we  arc  left  with  sound  speeds  that  arc  identical 
to  those  in  the  matrix  material.  Although  not  cited  by 
Carroll  and  Holt,  this  last  result  lends  more  weight  to  the 
arguments  in  favor  of  utilization  of  the  effective  stress 
model  ( Hq . (A.1)). 


APPENDIX  B 


11.  E.  TEST  PARAMETER  CALCULATIONS 


In  response  to  a request  by  the  contract  monitor, 

C.  R.  McFarland  of  Hdq. , DNA,  calculations  were  made  to 
investigate  the  relevance  of  material  properties  to  the 
ground  motion  and  stress  pulse  generated  by  high  explosive 
tests  in  NTS  tuff.  The  material  models  employed  were  based 
on  fits  to  laboratory  data  on  tuff  samples  obtained  from 
specific  locations  at  the  Nevada  Underground  Test  site. 

The  calculations  also  examined  the  sensitivity  of  the  results 
to  the  representation  of  the  high  explosive  energy  source  in 
the  calculations. 

Tuff  properties  from  two  locations  near  the  (U12e-12) 
underground  nuclear  test  site  were  considered  in  the  calcu- 
lations. Samples  of  the  tuff  at  distances  from  the  working 
point  of  30  ft  (Slifer  Hole  #1)  and  1330  ft  (Drill  Back  7) 

were  selected  from  those  tested  under  quasi-static  loading 

17  21  b 

by  Green,  et  al . These  two  locations  represent  the  ex- 

tremes in  the  measured  gas-filled  porosities  measured  at  the 
site.  At  the  Drill  BacK  7 (DB-7)  location  the  tuff  is 
almost  completely  saturated  with  only  1.6  percent  of  the 
volume  gas-filled.  At  Slifer  Hole  #1  (SH-1)  fully  7.6  per- 
cent of  the  volume  is  gas-filled,  representing  22%  of  the 
total  pore  volume  in  the  tuff.  The  solid-water-void  volume 
fractions  for  the  two  tuffs  are  listed  in  Table  B.l. 

A fit  to  the  quasi-static  data  from  the  DB-7  location 
was  used  by  Bjork^^  in  a predictive  calculation  for  a 
1000-lb  nitromethane  sphere  detonated  at  that  site.  The 
ground  motion  measurements  subsequently  made  in  thetest  were 
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TABLE  B.l 


DliSCRI PTION  OF  TEST  SPECIMENS 


SII-1 


l)B  - 7 


In -Situ  Density  (g/cc) 


1.86  + 0.03 


1.88  + 0.02 


luff  drain  Density  (g/cc)  2.35  + 


2.3  7 + 


v . 

drain  Volume  Fraction,  n 


0.05  5 


0.034 


’.Vatcr  Volume  Fraction,  n 


0.200 


0.350 


v > 

\ir  Volume  Fraction,  n 
(3)  / (2)  (3)  \ 


jl-J  l a J \ 
n / n + n 
0 \ 0 0 / 


0.070 


0.010 


4 . 4 °C 


in  good  agreement  with  the  calculations.  It  was  decided  to 
use  the  same  tuff  model  in  the  parameter  calculations  re- 
ported here  except  a treatment  of  the  irreversible  crushup 
states  was  added  in  order  to  permit  the  calculations  to 
follow  the  stress  pulse  propagation  to  greater  distances 
from  the  source  while  using  the  SKIPPFR  code. 

The  constitutive  relations  for  tuff  used  in  the  com- 
putations therefore  consisted  of  three  parts.  First,  the 
energy- pressure- volume  equation  of  state  for  the  completely 
crushed  material  is  specified  by  the  equation 


P = CpE  + Ap  + Bp2,  u = — 1 


where  the  reference  density  p here  refers  to  the  completely 

o 

crushed  tuff/water  mix  (only  water  filled  pores  remaining) 
at  ambient  conditions.  Second,  the  partially  crushed  material 
is  characterized  by  a porosity  parameter  ot (p ) = p (p)/p  (p)  . 

r 

Here  ps  and  p^  are  the  densities  of  the  completely  crushed 
and  partially  crushed  materials  respectively  at  the  same 
pressure.  A simple  form  for  a(p)  was  used  since  the  limited 
data  available  did  not  justify  a more  elaborate  fit, 

a = 1 + (a  -l)(l-p/pj2  (B.2) 

0 ^ 

Here  a is  the  value  of  a at  ambient  conditions  and  p 

o 

is  the  pressure  at  which  all  gas- filled  pores  are  crushed  out 
of  the  tuff.  The  p value  was  estimated  from  uniaxial 
strain  test  data. Third,  the  deviatoric  stress  was  com- 
puted with  an  elastic-perfectly  plastic  model,  using  a con- 
stant rigidity  modulus  G and  a simple  von  Mises  yield 
cond i t ion , 


S2  + S2  + S2  < 1 Y2  , (B.5) 

1 2 3 — •-> 

where  the  are  the  principal  deviatoric  stresses. 

The  values  of  the  constants  used  in  the  constitutive 
relations  for  the  two  partially  saturated  wet  tuffs  are 
listed  in  Table  B.2.  The  higher  volume  of  gas-filled  poro- 
sity at  the  SIl-1  location  may  be  expected  to  enhance  the 
displacement  at  given  pressure,  but  its  higher  shear  strength 
should  counteract  this  to  some  extent.  It  was  of  interest 
to  evaluate  these  competing  effects,  to  assist  in  the  selec- 
tion of  the  site  for  a planned  h.gh  explosive  test  in  the 
U16a  tunnel  complex  at  the  Nevada  Test  Site. 


w 


TABLE  B.2 

CONSTITUTIVE  MODEL  CONSTANTS  FOR 
TWO  LOCATIONS 


SH-1 

DB-7 

A (kbar) 

8 3.3 

83.3 

B (kbar) 

8 3.3 

83.3 

G 

1.8 

1.8 

P (g/cc) 
0 

2.012 

1.911 

a 

0 

1.082 

1.0165 

pc  (kbar) 

1.0 

0.5 

u (kbar) 

10.5 

10.5 

Y (kbar) 

0.55 

0.346 

The  source  was  taken  to  be  a 1000-lb  sphere  of 
Composition  B high  explosive  with  a density  of  =1.7  g/cc 

prior  to  detonation.  The  initial  radius  of  the  charge  is 

therefore  R = 39.5  cm.  Four  tuff  comparison  calculations 

o 

were  made  for  two  representations  of  the  source.  Two  cal- 
culations treated  the  burning  process  and  employed  a pressure- 
volume  -energy  equation  of  state  to  compute  the  detailed  be- 
havior of  the  detonation  products.  Two  additional  calculations 
were  made  in  which  the  source  was  approximated  by  a y-law 
gas  in  the  expanding  spherical  cavity  of  initial  radius  R . 

In  the  calculations  which  included  the  burning  pro- 
cess in  the  source  representation,  the  empirical  Jones- 

r 7 3 1 

Wilkins-Lee  equation  of  state1  1 was  used  to  describe  the 
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Composition  B detonation  products.  This  equation  is  pre- 
scribed by  the  pVF  relation 


The  constants  for  Composition  B are  listed  in  Table  B3.  In 
the  computations  the  detonation  front  propagates  from  the 
center  of  the  sphere  at  wave  speed  1)CJ  = 7.98  * 10 5 cm/sec  and 

the  chemical  energy  released  on  detonation  is  li  = 4.95  * 10 10 

o 

ergs/g.  The  corresponding  Chapman-Jouguet  pressure  and 
density  are  p^j  = 295  kbar  and  p^.j  = 2.35  g/cc,  respectively. 

In  the  two  calculations  in  which  the  source  was 
approximated  by  a cavity  containing  a y-law  gas,  the  pressure 
in  the  expanding  cavity  was  computed  from 

Y 

p ■ tB-5) 

where  the  value  y - 2.77  and  the  associated  parameters 
listed  in  Table  B.3  were  estimated  from  work  by  Coleburn 
and  Liddiard.  J The  density  and  pressure  are  assumed  to 
be  uniform  within  the  cavity  at  any  time.  The  initial 
cavity  pressure  in  the  computatir  s is  given  by 

2 7 7 

PQ  = 28^7T31t)  = 116  kbar  (B . 6) 

In  Fig.  B.l  the  profound  effect  of  air-filled  poro- 
sity on  shock  wave  attenuation  is  illustrated.  The  less 
saturated  SII-1  material  requires  only  about  half  the 
distance  to  attenuate  the  shock  to  a given  pressure  as 
does  the  DB-7  material.  The  y-law  gas  approximation  is 
seen  to  produce  results  very  close  to  those  obtained  using 
the  detailed  treatment  of  the  burning  process  in  the  source 
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tance  from  1000  lb  Composition 


TABLE  B . 3 


EQUATION  OF  STATE  PARAMETERS  FOR 
COMPOSITION  B HIGH  EXPLOSIVE: 


JWL 

y - Law 

A (dynes /cm2) 

5.24  x io  » 2 

- 

B (dynes/cm2) 

7.678  x io10 

- 

Cj  (g/cc) 

7.  21 

- 

C2  (g/cc) 

1.80 

- 

u> 

0.34 

- 

P (g/cc) 

1.717 

1.7 

Eo  (ergs/g) 

4.95  x io10 

- 

PCJ  (dynes/cm2) 

2.95  x io1 1 

2.83  x io11 

PCJ  Cg/cc) 

2.35 

2.35 

DCJ  (cm/sec)  j 

7.98  x 105 

7.95  x io 5 

representation.  At  early  times  the  y-law  approximation 
overestimates  the  driving  pressure  and  at  late  times  it 
underestimates  it.  This  is  illustrated  by  the  results  shown 
in  Fig . B . 2 . 

The  stress-time  history  at  radial  distances  from  the 
source  of  = 123.7  cm  and  203.6  cm,  respectively,  are 
shown  in  Figs.  B.3  and  B.4.  The  hump  in  the  profiles  for  the 
calculations  which  treated  the  detonation  processes  result 
from  a wave  which  is  reflected  at  the  center  of  the  cavity 
end  eventually  catches  up  with  the  shock  front. 

In  Fig.  B . 5 the  radial  displacement-time  history  is 
shown  at  a distance  of  R = 203.6  cm  from  the  source 
(Rq  - 207  cm  for  SH-1).  Although  the  calculations  have  not 


Burn-}  (a.  - 1.0165) 


Fig.  B4--Radial  stress-time  history  at  a radial  distance  ol 

R„  = 205.6  cm  from  the  1000  lb  Composition  15  hign  explosive 
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t 


been  carried  to  late  enough  times  to  determine  final  displace- 
ments, it  is  clear  that  much  greater  displacements  occur  in 
the  more  saturated  HR"  tuff.  At  this  distance  it  experiences 
a [leak  stress  of  5 whereas  the  SII-1  tuff  is  subjected 

to  a [leak  stress  of  only  ^ 1 kbar.  Finally,  the  time  of 
arrival  of  the  shock  front  at  a given  radial  distance  is 
shown  in  Fig.  R.6  for  each  of  the  two  tuffs. 
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Fig.  B6--Time  of  arrival  of  shock  front  at  indicated  radial  distance 
From  the  1000  lb  Composition  B high  explosive  source. 


APPliNDlX  C 
INTERACTION  TERMS 


(a)  In  the  Theor^  of  Interacting  Continue  (TINC) , each 
i has  a velocity  field  vJ  (x,t),  stress  field  (0a.  . rx  t) 

density  field  Vfx.t).  Associated  with  these^ixture 
quantities,  we  postulate  the  existence  of  (v,e  (x,t), 

°ij  (x,t)  and  p C fx>t)  in  the  physical  configuration. 

It  is  to  be  noted  that  neither  of  the  two  sets  of  quantities 
corresponds  to  the  actual  distribution  in  the  body.  The  quan- 
titles  in  the  physical  configuration  ( v)e,  ('p)e) 

represent  the  averages  of  the  actual  variables  tlken  over 
several  grains.  Thus,  purely  local  effects  (e.g.,  stress 
concentrations  along  a pore  boundary)  are  ignored.  The  re- 
lationships between  and  (gj®,  and  (g3  and  (p5e 

yej-e  given  in  Section  IV. 
v 


VJLJ  The  velocity  distributions 

(x.t)  and  v e (x.t)  will  not  be  in  general  identica 

In  the  following  discussion,  it  will  be,  however,  presumed 
that 


(a)  Cot) 

I)  _ D e 

0F~  = Ut~ 


whe  re 


(a) 

I) 

Dt 


(a) 

v • grad ; 


(a) 

v e *grad 


(C.l) 


This  follows  from  the  fact  that  one  is  trying  to  follow  the 
same  particle  in  both  the  mixture  and  the  physical  configuration. 

'e  following  discussion  will  be  restricted  to  a two- 
component  mixture.  Superscripts  (1)  and  (2)  are  used  to  de- 
note 5 (rock)  and  P (fluid),  respectively.  The  rock 
is  regarded  to  be  an  elastic-plastic  solid.  The  fluid  is 
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taken  to  be  inviscid,  but  due  account  is  taken  of  the  viscous 
effects  by  postulating  the  existence  of  a diffusive  force, 

-P  n.  It  represents  the  drag  force  experienced  by  the  fluid 

o — 

as  it  passes  through  the  porous  rock. 

(a) 

The  mass  conservation  relations  for  5 are: 


Mixture : 


(a) (a) 
n o 
fit — 


(a)  (a) 

+ p d i v v 


= 0 


ex  = 1,  2 


(C.2) 


Physical  Configuration: 


(‘VMc  (a)  (a) 

JLp^-  ♦ P d,,  v 


= 0 


a = 1,  2 


(C.  3) 


Combining  Eqs.  (1)  through  (3)  and  noting  the 


(a)  (a)  (a) 


n p 


(C.4) 


there  folljws: 


(a) 

div  V 


(a)  i 
div  V * jyy 

n 


r (a) 

9 n 

at 


(a) 

+ v 


( a) 

grad  n 


a = 1 , 2 


(C.5) 
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The  momentum  balance  relation  for 
is  given  by: 


(a) 

s in  the  mixture 


(a) 

P 


(a)  (a) 

D v. 
l 

Dt 


(a) 

9 Oh  Ca) 

“STp  * 0 Si 


a « 1,  2 


(c.(>) 


For  a two -component  mixture,  we  have 


CD  (2) 

PS  = -p  s = PS 


(C.  7) 


(2) 


For  6 , the  stress  tensor  a. . ( 

i.e.,  1J 


or 


Ca) 


a.?  ) 

ij  / 


ij  J 1S  isotropic, 


C2) 


(2) 


C2) (2) 


ij 


= * P «Aj  - * n P e 6,, 


C 2 ) (2) 


- J 


n o 


ij 


(C.  S) 


Noting  that  the  fluid  is  subjected  a drag  force  -p  n as  it 
passes  through  the  rock,  the  momentum  balance  relation  for  the 
fluid  in  the  physical  configuration  may  be  written  as : 


C2)eC2) 

D e Vie 

— nr1- 


C 2 ) 


dXi 


p n. 

o i 


CC.D) 


We  now  require  that 


'^eC2)e 
D v . 

I5t  1 


(2) (2) 
D v. 

l 

~Tt~ 


(C. 10) 


Relation  CC. 10)  follows  from  the  fact  that  the  change  of  momcn 
un  in  the  physical  configuration  must  equal  to  the  change  in 
the  mixture  configuration.  Therefore,  combining  Eqs . ((’  (,) 
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through  (10)  , we  obtain 


pc 


R = 


(2)c  (2) 

p grad  n 


(C.ll) 


Pi' 


consists  of  two  parts, 


Thus  the  momentum  interaction  term,  _ 
i.c,  dilatation  (-  P 0 grad  n ) and  shear  (p.n)  contribution*. 


Let  us  now 


(55 


\ amine  the  momentum  balance  relations  for 


(2)  \ 

the  rock  matrix , V.  In  contradistinction  to  the  fluid  I s . 
the  relationship  between  the  mixture  and  the  physical  configura- 
tions is  not  readily  obvious.  The  interaction  force,  pP,  may 
induce  a complex  state  of  stress  in  the  solid.  As  an  example, 
let  us  consider  the  mixture  to  he  in  a state  of  um-aijjal 
strain.  In  this  equation,  the  momentum  balance  for  s in 

the  mixture  is 


(1)(1) 

(1)  n v 

p 


U) 

3 o 


1 n 


i _ 


1 1 


»x 


(1) 

+ p R 


(C. 12) 


In  the  physical  configuration,  the  momentum  balance 
relation  fusing  (L.lOj)  is: 


(1) 


CD  CD  CD, 


n 


a a 


1 1 


Lit 


dX 


U), 

3 o 

i_ 

3x 


CD 

3 o 


i 3 


3x 


fC.  15) 


The  stress  components  Co’°  and  ll\'  will  in  general  be 
non-tero.  Although  the  mixture  is  undergoing  um-axinl  motion, 
the  individual  constituents  may  be  in  a more  (<jijmplex  stress 
state.  In  this  particular  case,  and  0)3  ® 

of  the  interactive  force,  pV.  (Combining  equations  (C.l*) 

( t . 13)  and  noting  that 


(a) 

a 


(a)  (a), 


(C.14) 


ij 


ij 
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there  follows: 


(1) 

p 8 = p8 

1 i 


CD 


= - a 


CD 

e 3 n 

T5T 


CD, 


3 a 


1 2 


1 1 


3x 


(u, 

3 a 
3x 


0 


From  Lqs.  (C.ll)  and  (C.15),  we  iiave  also: 


, . sT  t »>e 

ps,  ' ' p TT*  V,  = ’ °,1  — 


1 


1 


CD 

3 a 


CD, 

3 a 


1 2 


1 3 


3x 


3x 


0 


Thus,  in  trying  to  relate  th^^iixture  configuration  to  the 
physical  configuration  for  6 , it  is  necessary  t^^eparate 
out  the  stresses  induced  by  the  interaction  with  6 . Such 
an  identification  may  be  impossible  except  in  very  simple 
cases  such  as  uni-axial  motion. 


(a) 


The  internal  energy  balance  relation  for  i 
mixture  is 


in  the 


(a)  (a)  (a) 

(ct)  D E (a)  3 (a) 

0 Et  = °ij  "5*7  * p * 

a = 1 , 2 . 

We  will  now  assert  that: 

(a)e  (a) 

E e = E a = 1,  2 

(a) 

In  other  words,  the  internal  energy  E is  the  same  in  both 
the  mixture  and  the  physical  configurations. 


15) 


.;.io) 


:•  ID 


.18) 
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For 

tv  o parts  , 


the  fluid, 

(2) 

p 4 , and 


KCc« 

P 4' 


(2) 

11  assume  that  p 4* 


consists  of 


(2)  (2)  (2) 

o 4;  - P 4'j  + P Ts  (C.l'JJ 

(2)  . . (2) 

where  o : d denotes  the  dilatation  contribution  and  p 4S , 

the  diffusive  contribution. 

We  now  consider  the  energy  balance  equation  for  the 

fluid  in  the  physical  configuration 


(2) 

D 


e 


(2) (2) 
n 1; 
f) t 


(O.  (2)  (2) 

- p div  v + o ibs 


(C.  20) 


Combining  Lqs.  (C.4) 
wc  oil  tain: 


(2) 

P vd  = 


(C.5),  (C.S)  and  C C . 17)  througn  (C.20J 


(2)  (2) 
v • grad  n 


(C.21J 


Thus 


(2)  (2) 

3 • V 


(2) 


(2) 

= - p 


e 3 


(2) 


3 1 


- D 


(2) 

v 


(2) 


P 4w 


(C. 22) 


(1) 


For  the  solid,  5 , the  energy  balance  considerations 

are  much  more  complex.  F.ven  for  the  uni-axial  case  considered 

in  connection  with  the  momentum  balance,  the  situation  is 

1 1 J e , 

quite  difficult.  In  the  physical  configuration,  V2  and 
v C need  not  be  zero  and  indeed  may  depend  upon  all  the 
three  space  dimensions.  Fine , therefore,  needs  to  introduce 
some  additional  assumptions  for  the  solid  energy  balance. 

The  requirement  that  the  energy  contribution  of  the 
internal  material  interaction  forces  be  zero  may  be  written 
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as 


/Cl)  (2)  , 

P?  • I V - V ) 


(1)  (2) 

+ P 4>  + p ip  = 0 


C C . 2 3) 


Substituting  from  Hqs  . (C.ll),  (CMS)  and  (C.22)  into  liq.  (C  231 
there  follows:  4 1 J * 


C2)e  (2) 

P grad  n + p n 

o — 


(1)  Cl)  (2) 

V + p t/i  - p 


(2) 

e 3 n 


(2) 


<Tt 


p n * v 
o — 


(2) 

+ P = 0 


or 


/n’  C2), 


C2)  Cl) 


3 1 


P • grad  n * v 


Cl) 


+ i CD  C2)v 

+ p V + PQn  • ( v - v 


C2) 

+ P '1>S  = 0 


C C . 24) 


The  internal  energy  cent j-ibution  o'P  can  be  split  up  into 
two  parts,  P til'  and  PV  . 


Cl)  Cl)  Cl) 

p 'P  = P + p i/i 


CC. 25) 


This  yields 


C2) 


sT  (2>e  (2)  (1) 

Tt~  + p grad  n • v 


Cl)  Cl)  (2) 

+ P + p t/is  + p tiis 


/Cl)  (2)v 
+ P0H  * \ v ' v ) - 


(C. 26) 
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fl)  (1)  (2) 

Ke  have  thus  one  equation  for  three  unknowns , «l'd  . ^s*  's* 

We  will  now  assume  that  the  dilatation  and  shear  contributions 

must  he  separately  equal  to  tcro,  i.e., 


(1) 

P 


(2) 


(2) 


(2)  fl) 


1 ..  P—  + grad  n • v 


(C. 27a) 


(1)  (2) 


P + P = pnl  * ( V 


(2)  (1) 


(C. 27b) 


(1)  (2) 


I V x j v q\ 

This  procedure  is  exact  for  inyiscid  ( c ^ f 'Js  Po- 
or  incompressible  materials  ( n 3 constant).  However,  in 
the  general  case,  it  has  to  he  regarded  as  a constitutive 

assumption . 

The  next  question  that  arises  is  as  to  how  we  should 
split  the  diffusive  contribution.  Since  diffusion  is  a 
dissipative  process,  it  is  reasonable  to  require  that 


(1) 

P ’i's  1 & 


(2)  „ 
P 1 0 


One  further  assumption  is  now  necessary  regarding  the  par- 
tition of  the  diffusive  energy  contribution.  The  simplest 
hypothesis  is  to  assume  pVs  to  be  cero.  This  yields: 

(2)  / (2)  (l)v 

p «!»  = P n * v - v ) i 0 

s 0 ' 

The  last  assumption  merely  states  that  the  fluid  receiv 
all  the  diffusive  energy  contribution. 

In  our  previous  work  (3SR-648) , it  was  assumed  that 

(1)  t1) 

o lb  = -o  n * v 
rs  0—  - 


(C. 28) 


( C . 29) 


(2) 

p <l>  - p n 


(2) 

v 


(C. 30) 
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. 


1 


It  is  now  seen  that  such  an  assumption  cannot  possibly  satis- 
fy Eqs.  (C. 28) , and  hence  is  erroneous.  This  last  conclusion 
of  course  does  not  apply  if  p^n  is  regarded  as  a shear 
(non-dissipative)  force.  It  may  also  be  verified  that  for 
the  steady  case,  present  expressions  for  p ^ and  p ^ 
reduce  to  those  derived  in  3SR-648. 
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APPENDIX  D 

EFFECTIVE  STRESS  LAWS 


The  effect  of  fluid  pressure  on  the  deformation  and 
strength  of  saturated  rocks  has  been  studied  by  a number  of 
investigators.  This  role  of  fluid  pressure  is  usually  con- 


sidered by  defining  an  "effective  stress" 


where 


0.1) 


= applied  stress  on  the  rock/fluid  composite 
a = constant 


Pp  = pore  pressure 

The  'efiective  stress  law"  simply  implies  that  the  stress- 
strain  (and  strength)  response  of  the  saturated  rock  is 
identical  with  that  of  the  dry  rock  if  one  utilizes  the 
effective  stress  <aij>  instead  of  the  composite  stress 

°ij ' 

There  is  considerable  disagreement  as  to  the  value 

of  the  parameter  a.  Terzaghi^63^  argued  that  a should 

equal  the  porosity,  1-n,  so  that  the  effect  of  pore  pressure 

is  eliminated  when  the  porosity  is  zero.  Hubbert  and 

Rubey^  J,,61  attempted  to  prove  theoretically  that  a should 

be  1.  Although  the  validity  of  their  proof  is  somewhat 
fG  8 1 

controversial,1'  J experimental  measurements  on  strength 

have  generally  revealed  a fair  agreement  with  Eq . (D.l) 

17  7 -80  1 

with  a = 1. L J One  except  ion  to  this  good  agreement 

appears  to  be  in  the  area  of  strength  measurements  on  low 
porosity  rocks.  However,  Brace  and  Martin  ^80  ^ present 


Preceding  page  blank 
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quite  plausible  arguments  to  the  effect  that  this  apparent 
disagreement  is  due  to  fluid  pressure  in  not  being  fully 
effective.  Prior  to  gross  fracture,  the  rocks  begin  to  dilate 
under  compressive  stresses.  Above  a certain  critical  strain 
rate,  fluid  pressure  inside  the  pores  is  considerably  less 
than  the  applied  pore  pressure. 


I ven  though  a - 1 gives  good  agreement  with  strength 
tests,  its  application  to  s t res s - st 1 ain  measurements  is 
questionable. t81'8-5’68]  Skempton[82]  and  Gee  rlsina  [ 81  ] have 


suggested 

that 

a « 1 - K/Ks 

(D.2) 

wiie  re  K 

! bulk  modulus  of  porous  rock 

K 

s 

= bulk  modulus  of  rock  grain. 

Skulje[83  ] 

proposed  the  form 

a = 1 - n K/Kg  , 

(n.  3) 

where  n 

denotes  the  volume  fraction  of  the 

rock  matrix. 

Because  of  lack  of  good  theoretical  basis  or  sufficient  data, 
expressions  like  (D. 2)  and  (D.3)  have  not  been  generally 
employed.  Recently,  Nur  and  3yerlee^68^  have  derived  Eq. 
(D.2)  from  certain  linear  elasticity  considerations  and  tried 
to  correlate  it  with  s tress -strain  data. 

.>Iur  and  Byerlee^8’8^  have  tested  IVeber  sandstone 
under  hydrostatic  pressure  both  with  no  pore  fluid  and  when 
the  pores  are  saturated.  Choice  of  IVeber  sandstone  mini- 
mizes weakening  of  the  rock  matrix  by  chemico-phys ical 
interaction  with  the  pore  fluid.  They  measured  the  con- 
fining pressure  pc,  pore  pressure  p , and  the  partial 
volumetric  strain  for  the  rock,  0 . ‘in  Fig.  I).  1 we  show 
the  s t ress -strain  points  obtained  by  using  the  conventional 
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Strain  0 (10'3) 

Fig.  D. l--Conventional  effective  stress  law. 


ettective  stress  law.  It  is  seen  that  all  the  saturated 
data  lies  on  one  side  of  the  dry  data.  Had  the  conven- 
tional etfective  stress  law  been  applicable,  all  the  saturated 
data  should  have  fallen  on  the  dry  curve.  The  conventional 
effective  stress  law  overcorrocts  for  the  pore  pressure  effect 
and  is  evidently  inadequate.  In  Tig.  i),2  we  plot  the  stress- 
_-> 1 1 a 1 n data  obtained  from  the  Nur-Bycrlee  effective  stress 
lau.  Ihe  dry  stress-strain  data  was  numerically  differentiated 
(linear)  to  yield  K as  a discrete  function  of  p . Least 
square  technique  was  then  employed  to  express  K as  a second 
degree  polynomial  in  p . Following  \iur  and  Byerlee  ^ 8 ^ , 
hs  was  taken  to  be  0.36  mbar.  Tor  the  wet  case,  K was 
evaluated  from  the  polynomial  function  by  replacing  pr  by 
pc  ' pp ' Comparing  Figs.  D.  1 and  I).  2,  it  is  readily  seen 
that  use  of  the  Mur -fiver lee  effective  stress  law  considerably 
i educes  tiie  data  scatter.  However,  even  in  this  case  all  the 
wet  data  lie  on  one  side  of  the  dry  data.  Clearly,  the  Mur- 
fiyerlee  effective  stress  law  undercorrects  for  the  pore 
pressure  effect. 


In  Fig.  D.3,  we  plot  the  TINC  effective  stress  p e 
versus  effective  strain  0 e and  n.  The  interaction  function 
n was  then  expressed  as  a second  degree  polynomial  function 
in  pc.  For  the  wet  rock,  n was  evaluated  from  this  func- 
tional relationship  with  pc  replaced  by  pc  - 1.31  p . 
effective  pressure  p e and  effective  strain  ^0  e were 
then  determined  from  Eqs.  (3)  and  (23).  In  this  case,  tne 
scatter  of  wet  data  around  the  dry  data  is  considerably  re- 
duced. Also,  the  wet  data  lies  evenly  on  both  sides  of  the 
dry  data.  Thus  the  TINC  model  provides  a good  model  for 
data  fitting.  This,  however,  does  not  constitute  a proof 
of  the  model.  To  provide  a proof  of  the  model,  one  would 
also  need  to  measure  the  porosity  (1-n)  during  the  test  and 
compare  it  with  theoretically  calculated  values.  It  is, 
however,  comforting  to  observe  that  the  dependence  of  n. 


on 


and 


pp  is  in  accord  with  the  physical  intuition 


wet 
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